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Abstract 


Event and apparent horizons are key diagnostics for the presence and properties of 
black holes. In this article I review numerical algorithms and codes for Ending event 
and apparent horizons in numerically-computed spacetimes, focusing on calculations 
done using the 3-1-1 ADM formalism. 

The event horizon of an asymptotically-flat spacetime is the boundary between 
those events from which a future-pointing null geodesic can reach future null infin¬ 
ity, and those events from which no such geodesic exists. The event horizon is a 
(continuous) null surface in spacetime. The event horizon is defined nonlocally in 
time: it’s a global property of the entire spacetime, and must be found in a separate 
post-processing phase after (part of) the spacetime has been numerically computed. 

There are 3 basic algorithms for finding event horizons, based respectively on 
integrating null geodesics forwards in time, integrating null geodesics backwards in 
time, and integrating null surfaces backwards in time. The last of these is generally 
the most efficient and accurate. 

In contrast to an event horizon, an apparent horizon is defined locally in time in a 
spacelike slice, and so can be (and usually is) found “on the fly” during the numerical 
computation of a spacetime. A marginally outer trapped surface (MOTS) in a slice 
is a smooth closed 2-surface whose future-pointing outgoing null geodesics have zero 
expansion 0. An apparent horizon is then defined as a MOTS not contained in 
any other MOTS. The MOTS condition is a nonlinear elliptic partial differential 
equation (PDE) for the surface shape, containing the ADM 3-metric, its spatial 
derivatives, and the extrinsic curvature as coefficients. Most “apparent-horizon” 
finders actually find MOTSs. 

There are a large number of apparent-horizon finding algorithms, with differ¬ 
ing trade-offs between speed, robustness, accuracy, and ease of programming. In 
axisymmetry, shooting algorithms work well and are fairly easy to program. In 
slices with no continuous symmetries, Nakamura et aids algorithm and elliptic-PDE 
algorithms are fast and accurate, but require good initial guesses to converge. In 
many cases, Schnetter’s “pretracking” algorithm can greatly improve an elliptic- 
PDE algorithm’s robustness. Flow algorithms are generally quite slow, but can be 
very robust in their convergence. 
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Chapter 1 
Introduction 


Compact objects - ones which may contain event horizons and/or apparent 
horizons - are a major focus of numerical relativity. The usual output of a numerical 
relativity simulation is some (approximate, discrete) representation of the spacetime 
geometry (the 4-metric and possibly its derivatives) and any matter helds, but not 
any explicit information about the existence, precise location, or other properties of 
any event/apparent horizons. To gain this information, we must explicitly find the 
horizons from the numerically-computed spacetime geometry. The subject of this 
review is numerical algorithms and codes for doing this, focusing on calculations done 
using the 3-1-1 ADM formalism (iniiiii).^ Baumgarte and Shapiro [23 section 6] 
have also recently reviewed event and apparent-horizon hnding algorithms. 

In this review I distinguish between a numerical algorithm (an abstract descrip¬ 
tion of a mathematical computation; also often known as a “method” or “scheme”), 
and a computer code (a “horizon hnder”, a specihc piece of computer software 
which implements a horizon-hnding algorithm or algorithms). My main focus is on 
the algorithms, but I also mention specihc codes where they are freely available to 
other researchers. 

In this review I have tried to cover all the major horizon-hnding algorithms and 
codes, and to accurately credit the earliest publication of important ideas. However, 
in a held as large and active as numerical relativity, it’s inevitable that I have 
overlooked and/or misdescribed some important research. I apologise to anyone 
whose work I’ve slighted, and I ask readers to help make this a truly “living” review 
by sending me corrections, updates, and/or pointers to additional work (either their 
own or others) which I should discuss in future revisions of this review. 

The general outline of this review is as follows: In the remainder of this chapter 
I dehne notation and terminology fsection fl.lj) . discuss how 2-surfaces should be 

^ There are many interesting uses of event and/or apparent horizons in gaining physical under¬ 
standing of numerically-computed spacetimes. However, a discussion of these applications would 
encompass much of strong-held numerical relativity, and would be be far beyond the scope of this 
review. 
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parameterized fsection fl.2jl . and outline some of the software-engineering issues that 
arise in modern numerical relativity codes fsection ll.dj) . I then discuss numerical 
algorithms and codes for finding event horizons (chapter 12)) and apparent horizons 
(chapter inj). Finally, in the appendices I briefly outline some of the excellent nu¬ 
merical algorithms/codes available for two standard problems in numerical analysis, 
the solution of a single nonlinear algebraic equation (appendix El) and the time 
integration of a system of ordinary differential equations (appendix El) • 

1.1 Notation and Terminology 

I generally follow the sign and notation conventions of Wald jl4hj . I assume that 
all spacetimes are globally hyperbolic, and for event-horizon hnding I further assume 
asymptotic flatness; in this latter context is future null inhnity. I use the Penrose 
abstract-index notation, with summation over all repeated indices. 4-indices abc 
range over all spacetime coordinates {a;“}, and 3-indices ijk range over the spatial 
coordinates {x*} in a spacelike slice t = constant. The spacetime coordinates are 
thus = (t, x^). 

Qab is the spacetime 4-metric, and g°'^ the inverse spacetime 4-metric; these are 
used to raise and lower 4-indices. FJ^j, are the 4-Christoffel symbols. is the Lie 
derivative along the 4-vector held 

Qij is the 3-metric dehned in a slice, and is the inverse 3-metric; these are used 
to raise and lower 3-indices. V* is the associated 3-covariant derivative operator, 
and FF are the 3-Christoffel symbols, a and are the 3-1-1 lapse function and shift 
vector respectively,^ so the spacetime line element is 

ds‘^ = gabdx°'dx^ ( 1 . 1 ) 

= —{a^ — (3i(5'‘)dt^ + 2(3idx'‘dt + gijdx''dx^ (1.2) 

As is common in 3-1-1 numerical relativity, I follow the sign convention of Misner, 
Thorne, and Wheeler nnn and York jl49j in dehning the extrinsic curvature of the 
slice as Kij = —^Cngij = —Vinj, where is the future-pointing unit normal 
to the slice. (In contrast, Wald |145j omits the minus signs from this dehnition.) 
K = Ki is the trace of the extrinsic curvature Kij. ruADM is the ADM mass of an 
(asymptotically hat) slice. 

Indices uvw range over generic angular coordinates {9, (p) on S'^ or on a horizon 
surface. Note that these coordinates are conceptually completely distinct from the 
3-dimensional spatial coordinates a:*. Depending on the context, may or may 

not have the usual polar-spherical topology. 

Indices ijk label angular grid points on S '^ or on a horizon surface. Notice that 
these are 2-dimensional indices: a single such index uniquely specihes an angular 

^See York mni for a general overview of the 3-1-1 formalism as it’s used in numerical relativity. 
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grid point. 5ij is the Kronecker delta on the space of these indices, or equivalently 
on surface grid points. 

I often write a differential operator as F = F(|/, duV, duvV] Qij, dkQij^ Kij). The 
notation means that F is a (generally nonlinear) algebraic function of the variable y 
and its 1st and 2nd angular derivatives, and that F also depends on the coefficients 
Qij, dkQij, and Kij at the apparent horizon position. 

There are 3 common types of spacetimes/slices where numerical event or appar¬ 
ent horizon hnding is of interest: spherically-symmetric spacetimes/slices, axisym- 
metric spacetimes/slices, and spacetimes/slices with no continuous spatial symme¬ 
tries (no spacelike Killing vectors). I refer to the latter as “fully generic” space¬ 
times/slices. 

In this review I use the abbreviations “ODE” for ordinary differential equation, 
“PDE” for partial differential equation, “CE surface” for constant-expansion surface, 
and “MOTS” for marginally outer trapped surface. Names in Small Capitals 
refer to horizon hnders and other computer software. 

When discussing iterative numerical algorithms, it’s often useful to use the con¬ 
cept of an algorithm’s “radius of convergence”: Suppose the solution space within 
which the algorithm is iterating is S. Then given some norm || ■ || on S, the algorithm’s 
radius of convergence about a solution s G S is dehned as the smallest r > 0 such 
that the algorithm will converge to the correct solution s for any initial guess g with 
||g — s|| < r. We only rarely know the exact radius of convergence of an algorithm, 
but practical experience often provides a rough estimate.^ 

1.2 2-Surface Parameterizations 

1.2.1 Level-Set—Function Parameterizations 

The most general way to parameterize a 2-surface in a slice is to dehne a scalar 
“level-set function” F on some neighborhood of the surface, with the surface itself 
then being dehned as the level set 

F = 0 on the surface (1-3) 

Assuming the surface to be orientable, it’s conventional to choose F so that 
F > 0 (F < 0) outside (inside) the surface. 

This parameterization is valid for any surface topology, including time-dependent 
topologies. The 2-surface itself can then be found by a standard isosurface-hnding 
algorithm such as the marching-cubes algorithm [Mj. (This algorithm is widely used 

^An algorithm’s actual “convergence region” (the set of all initial guesses for which the algorithm 
converges to the correct solution) may even be fractal in shape. For example, the Julia set is the 
convergence region of Newton’s method on a simple nonlinear algebraic equation. 
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in computer graphics, and is implemented in a number of widely-available software 
libraries.) 

1.2.2 Strahlkorper Parameterizations 

Most apparent-horizon finders, and many event-horizon finders, assume that each 
connected component of the apparent (event) horizon has topology. With the 
exception of toroidal event horizons (discussed in section ITTl) . this is generally a 
reasonable assumption. 

To parameterize an surface’s shape, it’s common to further assume that we 
are given (or can compute) some “local coordinate origin” point inside the surface 
such that the surface’s 3-coordinate shape relative to that point is a “Strahlkorper”, 
(literally “ray body”, or more commonly “star-shaped region”) dehned by Minkowski 
f 11231 p. 108]) as 

a region in n-D Euclidean space containing the origin and whose surface, 
as seen from the origin, exhibits only one point in any direction. 

The Strahlkorper assumption is a significant restriction on the horizon’s coordi¬ 
nate shape (and the choice of the local coordinate origin). For example, it rules out 
the coordinate shape and local coordinate origin illustrated in figure 11.11 a horizon 
with such a coordinate shape about the local coordinate origin couldn’t be found by 
any horizon finder which assumes a Strahlkorper surface parameterization. 

For event-horizon Ending, algorithms and codes are now available which allow 
an arbitrary horizon topology, with no Strahlkorper assumption (see the discussion 
in section for details). For apparent-horizon finding, the flow algorithms 

discussed in section 13.2.71 theoretically allow any surface shape, although many im¬ 
plementations still make the Strahlkorper assumption. Removing this assumption 
for other apparent-horizon finding algorithms might be a fruitful area for further 
research. 

Given the Strahlkorper assumption, the surface can be explicitly parameterized 
as 

r = h{9, 0) (1.4) 

where r is the Euclidean distance from the local coordinate origin to a surface point, 
{6, 0) are generic angular coordinates on the horizon surface (or equivalently on S"^), 
and the “horizon shape function” h : ^ S?’*' is a positive real-valued function on 

the domain of angular coordinates defining the surface shape. 

There are two common discretizations of this surface representation: 

Spectral representation 

Here we expand the horizon shape function h in an infinite series in some 
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(typically orthonormal) set of basis functions such as spherical harmonics 
or symmetric trace-free tensors,"^ 

h{0, 0) = ^ aimYtmid, 4>) (1-5) 

£,m 

This series can then be truncated at some hnite order £max, and the iVcoefr = 
f'max(^max+l) Coefficients {aem} used to represent (discretely approximate) the 
horizon shape. For reasonable accuracy, f'max is typically on the order of 8 to 
12 . 

Finite difference representation 

Here we choose some hnite grid of angular coordinates {(6 *k,0k)}, k = 1, 2, 
3, ..., iVang on (or equivalently on the surface),® and represent (discretely 
approximate) the surface shape by the iVa ng values 

{h(6'K,0K)} K = 1, 2, 3, . . . , A^ang (1-6) 

For reasonable accuracy, iVa ng is typically on the order of a few thousand. 

^For convenience of exposition I use spherical harmonics here, but there are no essential differ¬ 
ences if other basis sets are used. 

discuss the choice of this angular grid in more detail in section 1 



Figure 1.1: This hgure shows a cross-section of a coordinate shape (the thick curve) 
which isn’t a Strahlkorper about the local coordinate origin shown (the large dot). 
The dashed line shows a ray from the local coordinate origin, which intersects the 
surface in more than one point. 


7 








It’s sometimes useful to explicitly construct a level-set function describing a given 
Strahlkorper. A common choice here is 

F = r-h{e,(p) (1.7) 

1.2.3 Finite-Element Parameterizations 

Another way to parameterize a 2-surface is via hnite elements, where the surface 
is modelled as a triangulated mesh, i.e. as a set of interlinked “vertices” (points 
in the slice, represented by their spatial coordinates {x*}), “edges” (represented 
by ordered pairs of vertices), and faces. Typically only triangular faces are used 
(represented as oriented triples of vertices). 

A key beneht of this representation is that it allows an arbitrary topology for the 
surface. However, determining the actual surface topology (e.g. testing for whether 
or not the surface self-intersects) is somewhat complicated. 

This representation is similar to that of Regge calculus mu ini and can 
similarly be expected to show 2nd order convergence with the surface resolution. 

Finite element surface representations have been used for apparent-horizon hnd- 
ing by Metzger 

1.3 Software-Engineering Issues 

Historically, numerical relativists wrote their own codes from scratch. As these 
became more complex, many researchers changed to working on “group codes” with 
multiple contributors. 

1.3.1 Software Libraries and Toolkits 

More recently, particularly in work on fully generic spacetimes, where all 3 spatial 
dimensions must be treated numerically, there has been a strong trend towards the 
use of higher-level software libraries and modular “computational toolkits” such as 
Cactus jHSI (http://www.cactuscode.org). These have a substantial learning 
overhead, but can allow researchers to work much more productively by focusing 
more on numerical relativity, instead of computer-science and software-engineering 
issues such as parameter-hle parsing, parallelization, I/O, etc. 

®There has been some controversy over whether, and if so how quickly, Regge calculus converges 
to the continuum Einstein equations. (See, for example, the debate between Brewin m and 
Miller mi, and the explicit convergence demonstration of Gentle and Miller EH]-) However, Brewin 
and Gentle m seem to have resolved this: Regge calculus does in fact converge to the continuum 
solution, and this convergence is generically 2nd order in the resolution. 



A particularly important area for such software infrastructure is mesh rehne- 
ment7 This is essential to much current numerical-relativity research, but is mod¬ 
erately difficult to implement even in only one spatial dimension, and much harder in 
multiple spatial dimensions. There are now a number of software libraries providing 
multi-dimensional mesh-rehnement infrastructure (sometimes combined with paral¬ 
lelization), such as DAGH/GrACE j1 07] (http://www.caip.rutgers.edu/~parashar/DAGH/) 
and ParaMesh (http: //ct .gsfc .nasa.gov/paramesh/Users_manual/ainr .html). 

The Cactus toolkit can be used in either unigrid or mesh-rehnement modes, the 
latter using a “mesh-rehnement driver” such as PAGH or Carpet [nni cm 
(http://www.carpetcode.org). 

In this review I point out event and apparent-horizon hnders which have been 
written in particular frameworks, and comment on whether they work with mesh 
rehnement. 

1.3.2 Code Reuse and Sharing 

Another issue is that of code reuse and sharing. It’s common for codes to be 
shared within a research group, but relatively uncommon for them to be shared 
between diherent (competing) research groups. Even apart from concerns about 
competitive advantage, without a modular structure and clear documentation it’s 
difficult to reuse another group’s code. The use of a common computational toolkit 
can greatly simplify such reuse. 

If such reuse can be accomplished, it becomes much easier for other researchers 
to build on existing work, rather than having to “reinvent the wheel”. As well as 
the obvious ease of reusing existing code that (hopefully!) already works and has 
been thoroughly debugged and tested, there’s another - less obvious - beneht of 
code sharing: It greatly eases the replication of past work, which is essential as a 
foundation for new development. That is, without access to another researcher’s 
code, it can be surprisingly difficult to replicate her results, because the success or 
failure of a numerical algorithm frequently depends on subtle implementation details 
not described in even the most complete of published papers. 

Event and apparent-horizon hnders are excellent candidates for software reuse: 

Many numerical-relativity researchers can beneht from using them, and they have 
a relatively simple interface to an underlying numerical-relativity simulation. Even 
if a standard computational toolkit isn’t used, this makes it relatively easy to port 
an event or apparent-horizon hnder to a diherent code. 

Throughout this review I note event and apparent-horizon hnders which are 
freely available to other researchers. 

^See, for example, Choptuik m, Pretorius and Pretorius and Choptuik uni for general 
surveys of the uses of, and methods for, mesh refinement in numerical relativity. 
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1.3.3 Using Multiple Event/Apparent Horizon Finders 

It’s often useful to have multiple event or apparent-horizon finders available: 
their strengths and weaknesses may complement each other, and the extent of 
(dis)agreement between their results can give a good measure of the numerical accu¬ 
racy. For example, figure EISl shows a comparison between the irreducible masses of 
apparent horizons in a binary black hole coalescence simulation (Alcubierre et al. 0 
hgure 4b]), as computed by two different apparent-horizon hnders in the Cactus 
toolkit, AHFinder and AHFinderDirect. In this case the two agree to within 
about 2% for the individual horizons, and 0.5% for the common horizon. 
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Chapter 2 

Finding Event Horizons 


2.1 Introduction 

The black hole region of an asymptotically-flat spacetime is defined (1721 EH) 
as the set of events from which no future-point null geodesic can reach future null 
inhnity {J^)- The event horizon is dehned as the boundary of the black hole region. 
The event horizon is a null surface in spacetime with (in the words of Hawking and 
Ellis [m p. 319]) “a number of nice properties” for studying the causal stucture of 
spacetime. 

The event horizon is a global property of an entire spacetime, and is dehned 
nonlocally in time: the event horizon in a slice is dehned in terms of (and can’t be 
computed without knowing) the full future development of that slice. 

In practice, to hnd an event horizon in a numerically-computed spacetime, we 
typically instrument a numerical evolution code to write out data hies of the 4- 
metric. After the evolution has reached an approximately-stationary hnal state, 
we then compute a numerical approximation to the event horizon in a separate 
post-processing pass, using the 4-metric data hies as inputs. 

As a null surface, the event horizon is necessarily continuous. In theory it need 
not be anywhere diherentiable,^ but in practice this behavior rarely occurs:^ The 
event horizon is generally smooth except for possibly a hnite set of “cusps” where 
new generators join the surface; the surface normal has a jump discontinuity across 
each cusp. (The classic example of such a cusp is the “inseam” of the “pair of pants” 
event horizon illustrated in hgures 12.31 and 12.41 ) 

A black hole is dehned as a connected component of the black hole region in 
a 3 -|- 1 slice. The boundary of a black hole (the event horizon) in a slice is a 

^Chrusciel and Galloway m showed that if a “cloud of sand” falls into a large black hole, each 
“sand grain” generates a non-differentiable caustic in the event horizon. 

^This is a statement about the types of spacetimes usually studied by numerical relativists, not 
a statement about the mathematical properties of the event horizon itself. 
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2-dimensional set of events. Usually this has 2-sphere (5^) topology. However, 
numerically simulating rotating dust collapse, Abrahams et al. P] found that in some 
cases the event horizon in a slice may be toroidal in topology. Lehner et al. [Sn] and 
Husa and Winicour |HI] have used null (characteristic) algorithms to give a general 
analysis of the event horizon’s topology in black hole collisions; they find that there 
is generically a (possibly brief) toroidal phase before the final 2-spherical state is 
reached. Lehner et al. ra later calculated movies showing this behavior for several 
asymmetric black hole collisions. 


2.2 Algorithms and Codes for Finding Event Hori¬ 
zons 

There are 3 basic event-horizon hnding algorithms: 

• Integrate null geodesics forwards in time fsection l2.2.1j) . 

• Integrate null geodesics backwards in time 1,section l2.2.2j) . 

• Integrate null surfaces backwards in time (section I2.2.3j) . 

I describe these in detail in the following subsections. 


2.2.1 Integrating Null Geodesics Forwards in Time 


The first generation of event-horizon finders were based directly on Hawking’s 
original definition of an event horizon: an event V is within the black hole region 
of spacetime if and only if there is no future-pointing “escape route” null geodesic 
from V to the event horizon is the boundary of the black hole region. 

That is, as described by Hughes et al. usi, we numerically integrate the null 
geodesic equation 


^ dx^ dx’^ 

*'(iA dX 


( 2 . 1 ) 


(where A is an affine parameter) forwards in time from a set of starting events and 
check which events have “escaping” geodesics. For analytical or semi-analytical 
studies like that of Bishop ra, this is an excellent algorithm. 

For numerical work it’s straightforward to rewrite the null geodesic equation m 
as a coupled system of two first-order equations, giving the time evolution of photon 
positions and 3-momenta in terms of the 3-1-1 geometry variables a, , and 

their spatial derivatives. These can then be time-integrated by standard numerical 
algorithms.^ However, in practice several factors complicate this algorithm: 

briefly review ODE integration algorithms and codes in appendix IBl 
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We typically only know the 3 + 1 geometry variables on a discrete lattice of 
spacetime grid points, and we only know the 3 + 1 geometry variables themselves, 
not their spatial derivatives. Therefore we must numerically differentiate the field 
variables, and interpolate the held variables and their spacetime derivatives to each 
integration point along each null geodesic. This is straightforward to implement,'^ 
but the numerical differentiation tends to amplify any numerical noise that may be 
present in the held variables. 

Another complicating factor is that the numerically computations generally only 
span a hnite region of spacetime, so it’s not entirely obvious whether or not a given 
geodesic will eventually reach . However, if the hnal numerically-generated slice 
contains an apparent horizon, we can use this as an approximation: any geodesic 
which is inside this apparent horizon will dehnitely not reach , while any other 
geodesic may be assumed to eventually reach if its momentum is directed away 
from the apparent horizon. If the hnal slice is approximately stationary, the error 
from this approximation should be small. (I discuss the “hnal slice is approximately 
stationary” assumption further in section 1^.2.3.11 1 

2.2.1.1 Spherically-Symmetric Spacetimes 

In spherical symmetry this algorithm works well, and has been used by a number 
of researchers. For example, Shapiro and Teukolsky used it to 

study event horizons in a variety of dynamical evolutions of spherically symmetric 
collapse systems. Figure shows an example of the event and apparent horizons 
in one of these simulations. 

2.2.1.2 Non—Spherically-Symmetric Spacetimes 

In a non-spherically-symmetric spacetime, several factors make this algorithm 
very inefficient: 

• Many trial events must be tried to accurately resolve the event horizon’s shape. 
(Hughes et al. jTH] describe a 2-stage adaptive numerical algorithm for choosing 
the trial events so as to accurately locate the event horizon as efficiently as 
possible.) 

• At each trial event we must try many different trial-geodesic starting directions 
to see if any of the geodesics escape to (or our numerical approximation to 
it). Hughes et al. |7H] report needing only 48 geodesics per trial event in several 
nonrotating axisymmetric spacetimes, but about 750 geodesics per trial event 
in rotating axisymmetric spacetimes, with up to 3000 geodesics per trial event 
in some regions of the spacetimes. 

'^In practice the differentiation can usefully be combined with the interpolation; I outline how 
this can be done in section Em 
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Figure 2.1: This figure shows part of a simulation of the spherically symmet¬ 
ric collapse of a model stellar core (a F = | polytrope) to a black hole. The 
event horizon (shown by the dashed line) was computed using the “integrate null 
geodesics forwards” algorithm described in section 12.2.11 solid lines show outgoing 
null geodesics. The apparent horizon (the boundary of the trapped region, shown 
shaded) was computed using the zero-finding algorithm discussed in section Id. 2. II 
The dotted lines show the world lines of Lagrangian matter tracers, and are labeled 
by the fraction of baryons interior to them. Figure reprinted with permission from 
Shapiro and Teukolsky, The Astrophysical Journal 235, 199-215 (1980), Copyright 
1980 by the American Astronomical Society. 
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• Finally, each individual geodesic integration requires many (short) time steps 
for an accurate integration, particularly in the strong-held region near the 
event horizon. 

Because of these limitations, for non-spherically-symmetric spacetimes the “in¬ 
tegrate null geodesics forwards” algorithm has generally been supplanted by the 
more efficient algorithms I describe in the following sections. 

2.2.2 Integrating Null Geodesics Backwards in Time 

It’s well-known that future-pointing outgoing null geodesics near the event hori¬ 
zon tend to diverge exponentially in time away from the event horizon. Figure EUl 
illustrates this behavior for Schwarzschild spacetime, but the behavior is actually 
quite generic. 

Anninos et al. |B] and Libson et al. [Oni observed that while this instability is 
a problem for the “integrate null geodesics forwards in time” algorithm (it forces 
that algorithm to take quite short time steps when integrating the geodesics), we 
can turn it to our advantage by integrating the geodesics backwards in time: the 
geodesics will now converge on to the horizon.^ 

This event-horizon finding algorithm is thus to integrate a large number of such 
(future-pointing outgoing) null geodesics backwards in time, starting on the hnal 
numerically-generated slice. As the backwards integration proceeds, even geodesics 
which started far from the event horizon will quickly converge to it. This can be 
seen, for example, in figures 12.11 and 12.21 

Unfortunately, this convergence property holds only for outgoing geodesics. In 
spherical symmetry the distinction between outgoing and ingoing geodesics is trivial, 
but as described by Libson et al. [22], 

[... ] for the general 3D case, when the two tangential directions of the 
EH are also considered, the situation becomes more complicated. Here 
normal and tangential are meant in the 3D spatial, not spacetime, sense. 
Whether or not a trajectory can eventually be “attracted” to the EH, and 
how long it takes for it to become “attracted,” depends on the photon’s 
starting direction of motion. We note that even for a photon which is 
already exactly on the EH at a certain instant, if its velocity at that 
point has some component tangential to the EH surface as generated by, 
say, numerical inaccuracy in integration, the photon will move outside of 
the EH when traced backward in time. For a small tangential velocity, 

®This convergence is only true in a global sense: locally the event horizon has no special geomet¬ 
ric properties, and the Riemann tensor components which govern geodesic convergence/divergence 
may have either sign. 
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Figure 2.2: This figure shows a number of light cones and future-pointing outgoing 
null geodesics in a neighborhood of the event horizon in Schwarzschild spacetime, 
plotted in ingoing Eddington-Finkelstein coordinates {t,r). (These coordinates are 
defined by the conditions that f -|- r is an ingoing null coordinate, while r is an areal 
radial coordinate.) The geodesics outside the event horizon are shown in green; those 
inside the event horizon are shown in red. All the geodesics start out close together 
near the event horizon; they diverge away from each other exponentially in time 
(here with an e-folding time of 4m near the horizon). Equivalently, they converge 
towards each other if integrated backwards in time (downwards on the page). 


16 





the photon will eventually return to the EH [... but] the position to 
which it returns will not be the original position. 

This kind of tangential drifting is undesirable not just because it in¬ 
troduces inaccuracy in the location of the EH, but more importantly, 
because it can lead to spurious dynamics of the “EH” thus found. Neigh¬ 
boring generators may cross, leading to numerically artihcial caustic 
points [...]. 


Libson et al. 


also observe that 


Another consequence of the second order nature of the geodesic equation 
is that not just the positions but also the directions must be specihed in 
starting the backward integration. Neighboring photons must have their 
starting direction well correlated in order to avoid tangential drifting 
across one another. 

Libson et al. 1331 give examples of the numerical difficulties that can result from 
these difficulties, and conclude that this event-horizon Ending algorithm 

[... ] is still quite demanding in finding an accurate history of the EH, 
although the difficulties are much milder than those arising from the 
instability of integrating forward in time. 


Because of this difficulty, this algorithm has generally been supplanted by the 
wards surface” algorithm I describe in the next section. 


‘back- 


2.2.3 Integrating Null Surfaces Backwards in Time 

Anninos et al. jHI, Libson et al. ra, and Walker mu introduced the important 
concept of explicitly (numerically) finding the event horizon as a null surface in 
spacetime. They observed that if we parameterize the event horizon with a (any) 
level-set function F satisfying (USD, then the condition for the surface F = 0 to be 
null is just 

g^^daFd,F = 0 ( 2 . 2 ) 

Applying a 3 -|- 1 decomposition to this then gives a quadratic equation which can 
be solved to find the time evolution of the level-set function. 


dtF = 


-g^^diF + ^J{g^^diFY - g^^g^^ifdjf 


9 


tt 


(2,3) 


Alternatively, assuming the event horizon in each slice to be a Strahlkorper in 
the manner of section imi we can define a suitable level-set function F by m- 
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Substituting this definition into (Q then gives an explicit evolution equation for 
the horizon shape function, 

o , _ -5'*'' + + \/(5'*'’ - g^'^duhy - - 2g^'^duh + g^^duhd^h) 

Surfaces near the event horizon share the same “attraction” property discussed in 
section 1^.2.21 for geodesics near the event horizon. Thus by integrating either surface 
representation (Q or (El backwards in time, we can refine an initial guess into a 
very accurate approximation to the event horizon. 

Notice that in contrast to the null geodesic equation EH), neither (El nor (El 
contain any derivatives of the 4-metric (or equivalently the 3 + 1 geometry variables). 
This makes it much easier to integrate these latter equations accurately.® 

This formulation of the event-horizon finding problem also completely eliminates 
the tangential-drifting problem discussed in section l2.2.2[ since the level-set function 
only parameterizes motion normal to the surface. 

2.2.3.1 Error Bounds: Integrating a Pair of Surfaces 

For a practical algorithm, it’s useful to integrate a pair of trial null surfaces 
backwards, an “inner-bound” one which starts (and thus always is) inside the event 
horizon, and an “outer-bound” one which starts (and thus always is) outside the 
event horizon. If the final slice contains an apparent horizon, then any 2-surface 
inside this can serve as our inner-bound surface. However, choosing an outer-bound 
surface is more difficult. 

It’s this desire for a reliable outer bound on the event horizon position that 
motivates our requirement for the final slice to be approximately stationary, since (in 
the absence of time-dependent equations of state or external perturbations entering 
the system) this ensures that, for example, any surface substantially outside the 
apparent horizon can serve as an outer-bound surface. 

Assuming we have an inner- and an outer-bound surface on the final slice, the 
spacing between these two surfaces after some period of backwards integration then 
gives a reliable error bound for the computed event horizon position. Equivalently, 
a necessary (and, if there are no other numerical problems, sufficient) condition for 
the event-horizon finding algorithm to be accurate is that the backwards integration 
must have proceeded far enough for the spacing between the two trial surfaces to be 
“small”. For a reasonable definition of “small”, this typically takes at least 15m adm 
of backwards integration, with 20m adm or more providing much higher accuracy. 

In some cases it’s difficult to obtain a long enough span of numerical data for 
this backwards integration. For example, in many simulations of binary black hole 

®Diener M describes how the algorithm can be enhanced to also determine (integrate) individ¬ 
ual null generators of the event horizon. This requires interpolating the 4-metric to the generator 
positions, but (still) not taking any derivatives of the 4-metric. 
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collisions, the evolution becomes unstable and crashes soon after a common apparent 
horizon forms. This means that we can’t compute an accurate event horizon for the 
most interesting region of the spacetime, that which is close to the black-hole merger. 
There’s no good solution to this problem except for the obvious one of developing a 
stable (or less-unstable) simulation that can be continued for a longer time. 

2.2.3.2 Explicit Strahlkorper Surface Representation 

The initial implementations of the “integrate null surface backwards” algorithm 
by Anninos et al. [H] , Libson et al. [12] , and Walker jl47j were based on the explicit 
Strahlkorper surface integration formula Q, further restricted to axisymmetry.^ 

For a single black hole the coordinate choice is straightforward. For the two- 
black-hole case, the authors used topologically cylindrical coordinates (p, z, 0), where 
the two black holes collide along the axisymmetry (z) axis. Based on the symmetry 
of the problem, they then assumed that the event horizon shape could be written in 
the form 

p = h{z) (2.5) 

in each t = constant slice. 

This spacetime’s event horizon has the now-classic “pair of pants” shape, with a 
non-differentiable cusp along the “inseam” (the z axis p = 0) where new generators 
join the surface. The authors tried two ways of treating this cusp numerically: 

• Since the cusp’s location is known a priori, it can be treated as a special 
case in the angular finite differencing, using one-sided numerical derivatives as 
necessary. 

• Alternatively, Thorne jl41j suggested calculating the union of the event hori¬ 
zon and all its null generators (including those which haven’t yet joined the 
surface). This “surface” has a complicated topology (it self-intersects along 
the cusp), but it’s smooth everywhere. This is illustrated by figure 1231 which 
shows a cross-section of this surface in a single slice, for a head-on binary black 
hole collision. For comparison, figure 12.41 shows a perspective view of part of 
the event horizon and some of its generators, for a similar head-on binary black 
hole collision. 

Caveny et al. 11011121 implemented the “integrate null surfaces backwards” algo¬ 
rithm for fully generic numerically-computed spacetimes, using the explicit Strahlkorper 
surface integration formula IQ. To handle moving black holes, they recentered 
each black hole’s Strahlkorper parameterization (d on the black hole’s coordinate 
centroid at each time step. 

^Walker mentions an implementation for fully generic slices, but only presents results for 
the axisymmetric case. 
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Figure 2.3: This figure shows a view of the numerically-computed event horizon in a 
single slice, together with the locus of the event horizon’s generators that haven’t yet 
joined the event horizon in this slice, for a head-on binary black hole collision. Notice 
how the event horizon is non-differentiable at the cusp where the new generators 
join it. Figure reprinted with permission from Libson et ai, Physical Review D 53, 
4335-4350 (1996). Copyright 1996 by the American Physical Society. 
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Figure 2.4: This figure shows a perspective view of the numerically-computed 
event horizon, together with some of its generators, for the head-on binary black 
hole collision discussed in detail by Matzner et al. Figure courtesy of Edward 
Seidel. 





For single-black-hole test cases (Kerr spacetime in various coordinates), they 
report typical accuracies of a few percent in determining the event horizon position 
and area. 

For binary-black-hole test cases (the Kastor-Traschen extremal-charge black hole 
coalescence with a cosmological constant), they detect black hole coalescence (which 
appears as a bifurcation in the backwards time integration) by the “necking off” of 
the surface. Figure EUl shows an example of their results. 

2.2.3.3 Level-Set Parameterization 

Caveny et al. EOl EH and Diener (independently) implemented the “in¬ 
tegrate null surfaces backwards” algorithm for fully generic numerically-computed 
spacetimes, using the level-set function integration formula dZSD. Here the level-set 
function F is initialized on the final slice of the evolution, and evolved backwards 
in time using (Q on (conceptually) the entire numerical grid. (In practice, only a 
smaller box containing the event horizon need be evolved.) 

This surface parameterization has the advantage that the event-horizon topology 
and (non-)smoothness are completely unconstrained, allowing the numerical study 
of conhgurations such as toroidal event horizons (discussed in section EE]). It’s also 
convenient that the level-set function F is defined on the same numerical grid as the 
spacetime geometry, so that no interpolation is needed for the evolution. 

The major problem with this algorithm is that during the backwards evolution, 
spatial gradients in F tend to steepen into a jump discontinuity at the event horizon,® 
eventually causing numerical difficulty. 

Caveny et al. |4()llTT] deal with this problem by adding an artificial viscosity term 
to the level-set function evolution equation, smoothing out the jump discontinuity 
in F. That is, instead of (ESD, they actually evolve F via 

dtF = e^V^F + RHS of Q (2.6) 

where is a generic 2nd order linear spatial differential operator, and £ > 0 is a 
(small) dissipation constant. This scheme works, but the numerical viscosity does 
seem to lead to significant errors (several percent) in their computed event-horizon 
positions and areas,® and even failure to converge to the correct solution for some 
test cases (e.g. rapidly-spinning Kerr black holes). 

Alternatively, Diener |33] developed a technique of periodically reinitializing the 
level-set function to approximately the signed distance from the event horizon. To 

^Equivalently, Diener observed that the locus of any given nonzero value of the level-set 
function F is itself a null surface, and tends to move (exponentially) closer and closer to the event 
horizon as the backwards evolution proceeds. 

®They describe how Richardson extrapolation can be used to improve the accuracy of the 
solutions from 0{e) to 0{e^), but it appears that this hasn’t been done for their published results. 
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Figure 2.5: This figure shows the cross-section of the numerically-computed event 
horizon in each of 5 different slices, for the head-on collision of two extremal Kastor- 
Traschen black holes. Figure reprinted with permission from Caveny and Matzner, 
Physical Review D 68, 104003 (2003). Copyright 2003 by the American Physical 
Society. 
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do this, he periodically evolves 


, ^ (|VF| - 1) 


(2.7) 


dxF = 


in an unphysical “pseudo-time” A until an approximate steady state has been achieved 
He reports that this works well in most circumstances, but can signihcantly distort 
the computed event horizon if the F = 0 isosurface (the current approximation to 
the event horizon) is only a few grid points thick in any direction, as typically occurs 
just around the time of a topology change in the isosurface. He avoids this problem 
by estimating the minimum thickness of this isosurface and, if it’s below a threshold, 
deferring the reinitialization. 

In various tests on analytical data, Diener EH found this event-horizon hnder, 
EHFinder, to be robust and highly accurate, typically locating the event horizon to 
much less than 1% of the 3-dimensional grid spacing. Even with only 10 grid points 
across the event horizon, this already corresponds to accuracies on the order of 0.1%, 
and this accuracy improves as expected (2nd order convergence) with increasing 
resolution. 

As an example of the results obtained with EHEinder, hgure 12.01 shows two 
views of the numerically-computed event horizon for a spiraling binary black hole 
collision. As another example, hgure IT7I shows the numerically-computed event and 
apparent horizons in the collapse of a rapidly rotating neutron star to a Kerr black 
hole. (The apparent horizons were computed using the AHFinderDirect code 
described in section umi) 

EHFinder is implemented as a module (“thorn”) in the Cactus computa¬ 
tional toolkit. It originally worked only with the PUGH unigrid driver, but Di¬ 
ener jSnj is currently enhancing it to work with the CARPET mesh-rehnement driver. 
EHFinder is freely available by anonymous CVS, and is now used by several 
research groups. 

2.3 Summary of Algorithms/Codes for Finding 
Event Horizons 

In spherical symmetry the “integrate null geodesics forwards” algorithm (sec¬ 
tion is reasonable, though the “integrate null geodesics backwards” algorithm 

fsection l2.2.2|l is more efficient. 

In non-spherically-symmetric spacetimes the “integrate null surfaces backwards” 
algorithm f.section l2.2.3j] is clearly the best algorithm known: it’s efficient, accurate, 
and fairly easy to implement. For generic spacetimes, Diener’s event-horizon hnder 
EHFinder is particularly notable as a freely available implementation of this 
algorithm as a module (“thorn”) in the widely-used Cactus computational toolkit. 
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Figure 2.6: This figure shows two views of the numerically-computed event horizon 
for a spiraling binary black hole collision. The initial data was constructed to have 
an approximate helical Killing vector, corresponding to black holes in approximately 
circular orbits (the D = 18 case of Grandclement et al. [TOj), with a proper separa¬ 
tion of the apparent horizons of 6.9m. Figure courtesy of Peter Diener, visualization 
by Werner Benger. 
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Figure 2.7: This figure shows the polar and equatorial radia of the event horizon 
(solid lines) and apparent horizon (dashed lines) in a numerical simulation of the 
collapse of a rapidly rotating neutron star to form a Kerr black hole. The dotted 
line shows the equatorial radius of the stellar surface. These results are from the D4 
simulation of Baiotti et al. ^H]- Notice how the event horizon grows from zero size, 
while the apparent horizon first appears at a finite size, and grows in a spacelike 
manner. Notice also that both surfaces are flattened due to the rotation. Figure 
reprinted with permission from Baiotti et al., Physical Review D 71, 024035 (2005), 


26 















Chapter 3 


Finding Apparent Horizons 

3.1 Introduction 

3.1.1 Definition 

Given a (spacelike) 3 + 1 slice, a “marginally outer trapped surface” (MOTS) is 
defined as a smooth (differentiable) closed orientable 2-surface in the slice whose 
future-pointing outgoing null geodesics have zero expansion 0. There may be 
multiple MOTSs in a slice. MOTSs may nest within each other, but they can’t 
cross. An “apparent horizon” is then dehned as an outermost MOTS in a slice, i.e., 
a MOTS not contained in any other MOTS. 

Equivalently, a “trapped surface” is dehned as a smooth closed 2-surface in the 
slice whose future-pointing outgoing null geodesics have negative expansion. The 
“trapped region” in the slice is then dehned as the union of all trapped surfaces, 
and the apparent horizon is dehned as the outer boundary of the trapped region. 

Notice that the apparent horizon is dehned locally in time (it can be computed 
using only Cauchy data on a spacelike slice), but (because of the requirement that it 
be closed) non-locally in space} Hawking and Ellis [71] discuss the general properties 
of MOTSs and apparent horizons in more detail. 

3.1.2 General Properties 

Given certain technical assumptions (including energy conditions), the existence 
of any MOTS (and hence any apparent horizon) implies that the slice contains 
a black hole.^ Moreover, the apparent horizon necessarily coincides with, or is 

^As an indication of the importance of the “closed” requirement, Hawking m has observed 
that if we consider two spacelike-separated events in Minkowski spacetime, the intersection of their 
backwards light cones satisfies all the conditions of the MOTS definition, except that it’s not closed. 

^Note that the converse of this latter statement is not true: An arbitrary (spacelike) slice 
through a black hole need not contain any apparent horizon. Notably, Wald and Iyer mni have 
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contained in, an event horizon. In a stationary spacetime the event and apparent 
horizons coincide. 

It’s this relation to the event horizon which makes apparent horizons valuable 
for numerical computation: an apparent horizon provides a useful approximation 
to the event horizon in a slice, but unlike the event horizon, an apparent horizon 
is dehned locally in time and so can be computed “on the fly” during a numerical 
evolution. 

Given a family of spacelike 3 + 1 slices which foliate part of spacetime, the 
union of the slices’ apparent horizons (assuming they exist) forms a world-tube. 
This world-tube is necessarily either null or spacelike. If it’s null this world-tube is 
slicing-independent (choosing a different family of slices gives the same world-tube, 
at least so long as each slice still intersects the world-tube in a surface with 2-sphere 
topology). However, if the world-tube is spacelike, it’s slicing-dependent: choosing 
a different family of slices will in general give a different world-tube.^ 

Except for flow algorithms (section I3.2.7|l , all numerical “apparent horizon” End¬ 
ing algorithms and codes actually find MOTSs, and hereinafter I generally follow the 
common (albeit sloppy) practice in numerical relativity of blurring the distinction 
between an MOTS and an apparent horizon. 

3.1.3 Trapping, Isolated, and Dynamical Horizons 

Hayward 1^ introduced the important concept of a “trapping horizon”, roughly 
speaking an apparent-horizon world-tube where the expansion becomes negative 
if the surface is deformed in the inward null direction, along with several useful 
variants. Ashtekar, Beetle, and Fairhurst jlSI and Ashtekar and Krishnan later 
defined the related concepts of an “isolated horizon”, essentially an apparent-horizon 
world-tube which is null, and a “dynamical horizon”, essentially an apparent-horizon 
world-tube which is spacelike. 

These world-tubes obey a variety of local and global conservation laws, and have 
many applications in analyzing numerically-computed spacetimes. See the references 
cited above, and also Dreyer et al. Ashtekar and Krishnan usiini, Gourgoulhon 
and Jaramillo cind Booth in2] for further discussions, including applications to 
numerical relativity. 

constructed a family of angularly anisotropic slices in Schwarzschild spacetime which approach 
arbitrarily close to r = 0 yet contain no apparent horizons. However, Schnetter and Krishnan ina 
have recently studied the behavior of apparent horizons in various anisotropic slices in Schwarz¬ 
schild and Vaidya spacetimes, finding that the Wald and Iyer behavior seems to be rare. 

^Ashtekar and Galloway m have recently proved “a number of physically interesting con¬ 
straints” on this slicing-dependence. 





3.1.4 Description in Terms of the 3 + 1 Variables 

In terms of the 3 + 1 variables, a marginally enter trapped surface (and thus an 
apparent horizon) satishes the condition f |148j . |721 section IIA]) 

Q = WiS^ + KijS^s^ -K = Q (3.1) 

where s* is the outward-pointing unit 3-vector normal to the surface.^ 

Assuming the Strahlkorper surface parameterization (El, (EH) can be rewritten 
in terms of angular 1st and 2nd derivatives of the horizon shape function h, 

O = 0(/l, Oy^hy dyyhy Qijy Kij'^ 0 (3.2) 

where 0 is a complicated nonlinear algebraic function of the arguments shown. 
(Shibata jl31j and Thornburg fTTninrij give the 0(h, dyh, dyyh) function explicitly.) 

3.1.5 Geometry Interpolation 

0 depends on the slice geometry variables gij, d^gij, and Kij at the horizon 
position.® In practice these variables are usually only known on the (3-dimensional) 
numerical grid of the underlying numerical-relativity simulation,® so they must be 
interpolated to the horizon position, and more generally, to the position of each 
intermediate-iterate trial shape the apparent-horizon hnding algorithm tries in the 
process of (hopefully) converging to the horizon position. 

Moreover, usually the underlying simulation gives only gij and Kij, so g^j must 
be numerically differentiated to obtain dkgij- As discussed by Thornburg moi sec¬ 
tion 6.1], it’s somewhat more efficient to combine the numerical differentiation and 
interpolation operations, essentially doing the differentiation inside the interpolator.^ 

Thornburg mni section 6.1] argues that for an elliptic-PDE algorithm (sec¬ 
tion 1X231), for best convergence of the nonlinear elliptic solver, the interpolated 

^Notice that in order for the 3-divergence in EU to be meaningful, s® (defined only as a field on 
the marginally outer trapped surface) must be smoothly continued off the surface, and extended 
to a field in some 3-dimensional neighborhood of the surface. The off-surface continuation is 
non-unique, but it’s easy to see that this doesn’t affect the value of 0 on the surface. 

®Or, in the Huq et al. 1231101 algorithm described in section n.2.5.21 at the local Cartesian grid 
point positions. 

®If the underlying simulation uses spectral methods (see Gottlieb and Orszag m and Boyd m 
for general discussions of spectral methods, and (for example) Ansorg et al. 13 El ^2, Bonazzola 
et al. [03E01I^, Grandclement et al. [03, Kidder et al. |13IHZIIH3, and Pfeiffer et al. |11()| for 
applications to numerical relativity) then the spectral series can be evaluated anywhere, so no 
actual interpolation need be done. 

^An interpolator generally works by (conceptually) locally fitting a fitting function (usually 
a low-degree polynomial) to the data points in a neighborhood of the interpolation point, then 
evaluating the fitting function at the interpolation point. By evaluating the derivative of the 
fitting function, the dkQij values can be obtained very cheaply at the same time as the gij values. 
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geometry variables should be smooth (differentiable) functions of the trial horizon 
surface position. He argues that that the usual Lagrange polynomial interpolation 
doesn’t suffice here (in some cases his Newton’s-method iteration failed to converge), 
because this interpolation gives results which are only piecewise differentiable.® To 
avoid this problem, Thornburg nm section 6.1] uses Hermite polynomial inter¬ 
polation; Cook and Abrahams M use bicubic spline interpolation. Most other 
researchers either don’t describe their interpolation scheme, or use Lagrange poly¬ 
nomial interpolation, and don’t report serious non-convergence problems. 

3.1.6 Criteria for Assessing Algorithms 

Ideally, an apparent-horizon finder should have several attributes: 

Robust: The algorithm/code should hnd an (the) apparent horizon in a wide range 
of numerically-computed slices, without requiring extensive tuning of initial 
guesses, iteration parameters, etc. This is often relatively easy to achieve for 
“tracking” the time evolution of an existing apparent horizon (where the most 
recent previously-found apparent horizon provides an excellent initial guess 
for the new apparent-horizon position), but may be difficult for detecting the 
appearance of a new (outermost) apparent horizon in an evolution, or for 
initial-data or other studies where there is no “previous time step”. 

Accurate: The algorithm/code should find an (the) apparent horizon to high ac¬ 
curacy and shouldn’t report spurious “solutions” (“solutions” which aren’t 
actually apparent horizons or, at least, marginally outer trapped surfaces). 

Efficient: The algorithm/code should be efficient in terms of its memory use and 
CPU time; in practice CPU time is generally the major constraint. It’s often 
desirable to find apparent horizons at each time step (or, at least, at frequent 
intervals) during a numerical evolution. For this to be practical the apparent- 
horizon finder must be very fast. 

In practice, no apparent-horizon finder is perfect in all these dimensions, so 
trade-offs are inevitable, particularly when ease of programming is considered. 

As discussed in section oi there are also significant advantages to having an 
apparent-horizon finder that’s freely available to other research groups, particularly 
if it’s designed and documented in such a way as to be relatively portable. 

3.1.7 Local versus Global Algorithms 

Apparent-horizon finding algorithms can usefully be divided into two broad 
classes: 

^Thornburg [HHl appendix F] gives a more detailed discussion of this non-smoothness of 
Lagrange-polynomial interpolation errors. 
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Local algorithms are those whose convergence is only gnaranteed in some (fnnc- 
tional) neighborhood of a solntion. These algorithms reqnire a “good” initial 
gness in order to find the apparent horizon. Most apparent-horizon finding 
algorithms are local. 

Global algorithms are those which can (in theory, ignoring finite-step-size and 
other nnmerical effects) converge to the apparent horizon independent of any 
initial guess. Flow algorithms (section ld.2.7|l are the only truely global al¬ 
gorithms. Zero-finding in spherical symmetry (section ld.2.1|l and shooting 
in axisymmetry (section ld.2.2|l are “almost global” algorithms: they require 
only 1-dimensional searches, which (as discussed in appendix can be pro¬ 
grammed to be very robust and efficient. In many cases horizon pretracking 
(section I3.2.fijl can semi-automatically find an initial guess for a local algo¬ 
rithm, essentially making the local algorithm behave like an “almost global” 
one. 

One might wonder why local algorithms are ever used, given the apparently 
superior robustness (guaranteed convergence independent of any initial guess) of 
global algorithms. There are two basic reasons: 

• In practice, local algorithms are much faster than global ones, particularly 
when “tracking” the time evolution of an existing apparent horizon. 

• Due to finite-step-size and other numerical effects, in practice even “global” 
algorithms may fail to converge to an apparent horizon (that is, the algorithms 
may sometimes fail to find an apparent horizon even when one exists in the 
slice). 


3.2 Algorithms and Codes for Finding Apparent 
Horizons 

Many researchers have studied the apparent-horizon-finding problem, and there 
are a large number of different apparent-horizon finding algorithms and codes. 
Almost all of these require (assume) that any apparent horizon to be found is a 
Strahlkorper (section ll.2|l about some local coordinate origin; both finite-difference 
and spectral parameterizations of the Strahlkorper are common. 

For slices with continuous symmetries, special algorithms are sometimes used: 

Zero-Finding in Spherical Symmetry fsection ld.2.1j) 

In spherical symmetry the apparent horizon equation (EH) becomes a 1-dimensional 
nonlinear algebraic equation, which can be solved by zero-finding. 
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The Shooting Algorithm in Axisymmetry fsection l;i.2.2jl 

In axisymmetry the apparent horizon equation (EH) becomes a nonlinear 2- 
point boundary value ODE, which can be solved by a shooting algorithm. 

Alternatively, all the algorithms described below for generic slices are also appli¬ 
cable to axisymmetric slices, and can take advantage of the axisymmetry to simplify 
the implementation and boost performance. 

For fully generic slices, there are several broad categories of apparent-horizon 
Ending algorithms and codes: 

Minimization Algorithms fsection Id. 2. dj) 

These algorithms dehne a scalar norm on 0 over the space of possible trial 
surfaces. A general-purpose scalar-function-minimization routine is then used 
to search trial-surface-shape space for a minimum of this norm (which should 
give 0 = 0). 

Nakamura et aZ.’s Spectral Integral-Iteration Algorithm fsection rd.2.4jl 

This algorithm expands the (Strahlkorper) apparent-horizon shape function 
in a spherical-harmonic basis, uses the orthogonality of spherical harmonics 
to write the apparent horizon equation as a set of integral equations for the 
spectral coefficients, and solves these equations using a functional-iteration 
algorithm. 

Elliptic-PDE Algorithms fsection Fd. 2. 5|1 

These algorithms write the apparent horizon equation (EH as a nonlinear 
elliptic (boundary-value) PDE for the horizon shape, and solve this PDF using 
(typically) standard elliptic-PDE numerical algorithms. 

Horizon Pretracking fsection ld.2.b|l 

Horizon pretracking solves a slightly more general problem than apparent- 
horizon finding: roughly speaking, the determination of the smallest E > 0 
such that the equation 0 = E has a solution, and the determination of that 
solution. By monitoring the time evolution of E and of the surfaces satisfying 
this condition, we can determine - before it appears - approximately where 
(in space) and when (in time) a new marginally outer trapped surface will 
appear in a dynamic numerically-evolving spacetime. Horizon pretracking 
is implemented as a 1-dimensional (binary) search using a slightly-modified 
elliptic-PDE apparent-horizon finding algorithm as a “subroutine”. 

Flow Algorithms (section EHZl) 

These algorithms start with a large 2-surface (larger than any possible ap¬ 
parent horizon in the slice), and shrink it inwards using an algorithm which 
ensures that the surface will stop shrinking when it coincides with the apparent 
horizon. 
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I describe the major algorithms and codes in these categories in detail in the 
following subsections. 

3.2.1 Zero-Finding in Spherical Symmetry 

In a spherically symmetric slice, any apparent horizon must also be spherically 
symmetric, so the apparent horizon equation (HI becomes a 1-dimensional nonlin¬ 
ear algebraic equation 0(h) = 0 for the horizon radius h. For example, assuming the 
usual polar-spherical spatial coordinates x* = (r, 9, 0), we have 1 |ld8| equation (B7)]) 

0 = ^ Q 

geey/^ gee 

Given the geometry variables grr, gee, dr gee, and Kee, this equation may be easily and 
accurately solved using one of the zero-hnding algorithms discussed in appendix 1X1^ 
Zero-hnding has been used by many researchers, including inni inzi inn 1123 
I1U91 imi Il24l IHl 1189] .^° For example, the apparent horizons shown in hEfure I^TD 
were obtained using this algorithm. As another example, hgnre lim shows 0(r) and 
h at various times in a (different) spherically symmetric collapse simulation. 

3.2.2 The Shooting Algorithm in Axisymmetry 

In an axisymmetric spacetime, the space of angular coordinates (6^, 0) is effec¬ 
tively 1-dimensional, and given the Strahlkorper assumption, without further loss of 
generality we can write the horizon shape function as h = h{6), where 6 is the single 
nontrivial angular coordinate. The apparent horizon equation (HU) then becomes 
a nonlinear 2-point boundary-value ODE for the horizon shape function h (USD 
equation (1.1)]) 

O = 0(/i, dehy deeh, gij, dkgtj, ^ij') 0 (^•^) 

where 0(h) is a nonlinear 2nd order (ordinary) differential operator in h as shown. 

Taking the angular coordinate 6 to have the usual polar-spherical topology, local 
smoothness of the apparent horizon gives the boundary conditions 

deh = 0 at 9=0 and 9=9 ^^^ (3.5) 

where 6*max is vr/2 if there is “bitant” reflection symmetry about the z = 0 plane, or 
71 otherwise. 

®Note that drQee is a known coefficient field here, not an unknown (if necessary, it can be 
obtained by numerically differentiating gge). Therefore, despite the appearance of the derivative, 
is still an algebraic equation for the horizon radius h, not a differential equation. 

^°See also the work of Bizon, Malec, and O Murchadha |SH1 for an interesting analytical study 
giving necessary and sufficient conditions for apparent horizons to form in non-vacuum spherically 
symmetric spacetimes. 
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Figure 3.1: This figure shows results for a spherically symmetric numerical evolution 
of a black hole accreting a narrow shell of scalar field, the SOO.pqwl evolution of 
Thornburg [na. Part (a) of this figure shows 0(r) (here labelled H) for a set 
of equally-spaced times between t=19 and t=20, while part (b) shows the corre¬ 
sponding horizon radius h{t) and the Misner-Sharp jlOOj mass m{h) internal to each 
marginally outer trapped surface (MOTS). Notice how two new MOTSs appear 
when the local minimum in 0(r) touches the 0=0 line, and two existing MOTS 
disappear when the local maximum in 0(r) touches the 0=0 line. 
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As well as the more general algorithms described in the following snbsections, 
this may be solved by a shooting algorithm: 

1. Gness the valne of h at one endpoint, say h{9=0) = h*. 

2. Use this gnessed valne of h{9=0) together with the bonndary condition there ()3.5|1 
as initial data to integrate (“shoot”) the ODE (jSill) from 9=0 to the other 
endpoint 6*=6'max-^^ 

3. If the nnmerically compnted solntion satisfies the other bonndary condition (unD 
at 9=9jnax to within some tolerance, then the jnst-computed h{9) describes the 
(an) apparent horizon, and the algorithm is finished. 

4. Otherwise, adjnst the gnessed valne h{9=0) = and try again. Because 
there’s only a single parameter (h*) to be adjusted, this can be done easily 
and efficiently using one of the 1-dimensional zero-finding algorithms discussed 
in appendix 1X1 

This algorithm is fairly efficient and easy to program. By trying a sufficiently 
wide range of initial guesses h* this algorithm can give a high degree of confidence 
that all apparent horizons have been located, although this, of course, increases the 
cost. 

Shooting algorithms of this type have been used by many researchers, for example 

[m Ea El Ea EEi cnni El n . 

3.2.3 Minimization Algorithms 

This class of algorithms defines a scalar norm || ■ || on the expansion 0 over the 
space of possible trial surfaces, typically the mean-squared norm 

02 dU (3.6) 

where the integral is over all solid angles on a trial surface. 

Assuming the horizon surface to be a Strahlkorper and adopting the spectral 
representation for the horizon surface, we can view the norm (EH) as being 

defined on the space of spectral coefficients {aim}- 

This norm clearly has a global minimum ||0|| = 0 for each solution of the appar¬ 
ent horizon equation dn. To find the apparent horizon we numerically search the 
spectral-coefficient space for this (a) minimum, using a general-purpose “function- 
minimization” algorithm (code) such as Powell’s algorithm. 

briefly review ODE integration algorithms and codes in appendix 0 
^^See, for example, Dennis and Schnabel m or Brent m for general surveys of general-purposes 
function-minimization algorithms and codes. 
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Evaluating the norm (ESI) requires a numerical integration over the horizon 
surface: We choose some grid of Wng points on the surface, interpolate the slice 
geometry helds {gij, dkQij, and Kij) to this grid (see section I!?. 1.5j) . and use numerical 
quadrature to approximate the integral. In practice this must be done for many 
different trial surface shapes (see section in.2.d.2jl . so it’s important that it be as 
efficient as possible. Anninos et al. jlj and Baumgarte et al. discuss various 
ways to optimize and/or parallelize this calculation. 

Unfortunately, minimization algorithms have two serious disadvantages for apparent- 
horizon Ending: they are susceptible to spurious local minima, and they’re very slow. 

I discuss these disadvantages further in the following two subsections. 

3.2.3.1 Spurious Local Minima 

While the norm dSH) clearly has a single global minimum || 0 || = 0 for each 
marginally outer trapped surface 0 = 0 , it typically also has a large number of 
other local minima with 0 7 ^ 0 , which are “spurious” in the sense that they don’t 
correspond (even approximately) to marginally outer trapped surfaces. Unfortu¬ 
nately, general-purpose “function-minimization” routines only locate local minima, 
and thus may easily converge to one of the spurious 0 7 ^ 0 minima. 

What this problem means in practice is that a minimization algorithm needs 
quite a good (accurate) initial guess for the horizon shape in order to ensure that 
the algorithm converges to the true global minimum 0 = 0 rather than to one of 
the spurious 0 7 ^ 0 local minima. 

To view this problem from a difierent perspective, once the function-minimization 
algorithm does converge, we must somehow determine whether the “solution” found 
is the true one 0 = 0 or a spurious one 0 7 ^ 0. Due to numerical errors in the 
geometry interpolation and the evaluation of the integral ()3.f)|l . ||0|| will almost 
never evaluate to exactly zero; rather, we must set a tolerance level for how large 
||0|| may be. Unfortunately, in practice it’s hard to choose this tolerance: if it’s 
too small, the genuine solution may be falsely rejected, while if it’s too large, we 
may accept a spurious solution (which may be very diEerent from any of the true 
solutions). 

Anninos et al. [7] and Baumgarte et al. [221 suggest screening out spurious 
solutions by repeating the algorithm with varying resolutions of the horizon-surface 
grid, and checking that ||0|| shows the proper convergence towards zero. This seems 
like a good strategy, but it’s tricky to automate and, again, it may be difficult to 

^^There’s a simple heuristic argument (see, for example, Press et al. mu section 9.6]) that at 
least some spurious local minima should be expected: We are trying to solve a system of Wng 
nonlinear equations ©i = 0 (one equation for each horizon-surface grid point). Equivalently, we are 
trying to find the intersection of the Wng codimension-one hypersurfaces ©i = 0 in surface-shape 
space. The problem is that anywhere two or more of these hypersurfaces closely approach, but 
don’t actually touch, there is a spurious local minimum in ||©||. 
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choose the necessary error tolerances in advance. 


3.2.3.2 Performance 

For convenience of exposition, snppose the spectral representation (EH) of the 
horizon-shape fnnction h nses spherical harmonics (Symmetric trace-free ten¬ 
sors or other basis sets don’t change the argnment in any important way.) Then if 
we keep harmonics np to some maximnm degree ^max, the nnmber of coefficients is 
then iVcoeff = (t'max+1)^- -^max IS Set by the desired accuracy (angular resolution) of 
the algorithm, and is typically on the order of 6 to 12. 

To hnd a minimum in an Atcogfr-dimensional space (here the space of surface- 
shape coefficients {aim}), a general-purpose function-minimization algorithm typ¬ 
ically needs on the order of to lOiV^^gg iterations. Thus the number of 

iterations grows as f'^ax- 

Each iteration requires an evaluation of the norm (EH for some trial set of 
surface-shape coefficients {aim}, which requires 0{Ncoes) = Q(^ ma x) work to com¬ 
pute the surface positions, together with O(iVang) work to interpolate the geometry 
fields to the surface points and compute the numerical quadrature of the inte¬ 
gral (EH- 

The result is that minimization horizon-hnders tend to be quite slow, particu¬ 
larly if high accuracy is required (large £max and A^ng)- The one exception is in 
axisymmetry, where only spherical harmonics Yim with m=0 need be considered. 
In this case minimization algorithms are much faster, though probably still slower 
than shooting or elliptic-PDE algorithms. 

simple counting argument suffices to show that any general-purpose function-minimization 
algorithm in n dimensions must involve at least 0{n^) function evaluations (see, for example, Press 
et al. irm section 10.6]): Suppose the function to be minimized is / : Jff" ^ 3?, and suppose / 
has a local minimum near some point xq G 3?". Taylor-expanding / in a neighborhood of xo gives 
/(x) = /(xo) -I- a^(x—Xo) -I- (x—xo)^B(x—Xo) -I- 0(||x—xojp), where a G 3ff”, B G 3?”’^” is symmetric, 
and denotes the transpose of the column vector v G 3?". 

Neglecting the higher order terms (i.e., approximating / as a quadratic form in x in a neighbor¬ 
hood of Xo), and ignoring /(xo) (which doesn’t affect the position of the minimum), there are a 
total oi N = n + ^n{n + 1) coefficients in this expression. Changing any of these coefficients may 
change the position of the minimum, and at each function evaluation the algorithm “learns” only 
a single number (the value of / at the selected evaluation point), so the algorithm must make at 
least N = 0{n‘^) function evaluations to (implicitly) determine all the coefficients. 

Actual functions aren’t exact quadratic forms, so in practice there are additional 0(1) multi¬ 
plicative factors in the number of function evaluations. Minimization algorithms may also make 
additional performance and/or space-versus-time trade-offs to improve numerical robustness or to 
avoid explicitly manipulating n x n Jacobian matrices. 
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3.2.3.3 Summary of Minimization Algorithms/Codes 

Minimization algorithms are fairly easy to program and have been used by many 
researchers, for example 123 123 123 13 E3 El- However, they’re susceptible to 
spurious local minima, have relatively poor accuracy, and tend to be quite slow. 
I believe that the other algorithms discussed in the following sections are generally 
preferable. 

Alcubierre’s apparent-horizon hnder AHFinder j3] includes a minimization al¬ 
gorithm based on the work of Anninos et al. |3-^^ It’s implemented as a module 
(“thorn”) in the CACTUS computational toolkit, and is freely available by anony¬ 
mous CVS (it’s part of the CactusEinstein set of thorns included with the stan¬ 
dard Cactus distribution). It has been used by a number of research groups. 

3.2.4 Nakamura et aL's Spectral Integral-Iteration Algo¬ 
rithm 

Nakamura, Kojima, and Oohara [102] developed a functional-iteration algorithm 
for solving the apparent horizon equation (HI. 

This algorithm begins by choosing the usual polar-spherical topology for the 
angular coordinates (6*, 0), and rewriting the apparent horizon equation ()3.2|1 in the 
form 

L = deeh + p\ + ^ = F{de^h, deh, d^h; TJ,-) (3.7) 

tcLii C7 sm u 

where F is a complicated nonlinear algebraic function of the arguments shown, which 
remains regular even at 6=0 and 6=7r, and where for future use we dehne L to be 
the left hand side of (II3)- 

Next we expand h in spherical harmonics Because the left hand side L 

of ()3.7|l is just the flat-space angular Laplacian of h, which has the as orthogonal 
eigenfunctions, multiplying both sides of (ISI3 by (the complex conjugate of Y^m) 
and integrating over all solid angles gives 

J (3-8) 

for each i and m except i = m = 0. 

Based on this, Nakamura et al. [102] propose the following functional-iteration 
algorithm for solving (ISI3: 

1. Start with some (initial-guess) set of horizon-shape coefficients {aim}- These 
determine a surface shape via (EH). 

2. Interpolate the geometry variables to this surface shape (see section IH.l.bj) . 
AHFinder also includes a “fast flow” algorithm ('section B.2.7II . 
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3. For each ^ and m except £ = m = 0, evaluate the integral (ESI) by numerical 
quadrature to obtain a next-iteration coefficient agm- 

4. Determine a next-iteration coefficient aoo by numerically solving (hnding a 
root of) the equation 

jY*QFdn = 0 (3.9) 

This can be done using any of the 1-dimensional zero-hnding algorithms dis¬ 
cussed in appendix 1X1 

5. Iterate until all the coefficients {agm} converge. 

Gundlach 112 ] observed that the subtraction and inversion of the flat-space 
angular Laplacian operator in this algorithm is actually a standard technique for 
solving nonlinear elliptic PDFs by spectral methods. I discuss this observation and 
its implications further in section 13.2.7.41 

Nakamura et al. report that their algorithm works well, but Nakao mu 
has argued that it tends to become inefficient (and possibly inaccurate) for large 
i (high angular resolution) because the Yg^ fail to be numerically orthogonal due 
to the hnite resolution of the numerical grid. I know of no other published work 
addressing Nakao’s criticism. 

Kemball and Bishop IMI investigated the behavior of Nakamura et aids algo¬ 
rithm, and found that its (only) major weakness seems to be that the Ooo-update 
equation (ESD “may have multiple roots or minima even in the presence of a single 
marginally outer trapped surface, and all should be tried for convergence”. 

Kemball and Bishop jSl] suggested and tested several modihcations to improve 
the algorithm’s convergence behavior. They verihed that (either in its original form 
or with their modihcations) the algorithm’s convergence speed (number of iterations 
to a given error level) is roughly independent of the degree fmax of spherical-harmonic 
expansion used. They also give an analysis that the algorithm’s cost is and 

its accuracy £ = 0(l/^max), be. the cost is 0{l/e^). 

Despite what appears to be fairly good numerical behavior and reasonable ease 
of implementation, this algorithm has not been widely used apart from later work 
by its original developers (see, for example, 

3.2.5 Elliptic-PDE Algorithms 

The basic concept of elliptic-PDE algorithms is simple: we view the apparent 
horizon equation (EH as a nonlinear elliptic PDF for the horizon shape function h on 
the angular-coordinate space and solve this equation by standard hnite-differencing 
techniques,^® generally using Newton’s method to solve the resulting set of nonlinear 

theory this equation could also be solved by a spectral method on using spectral 
differentiation to evaluate the angular derivatives. (See the references cited in footnoteElon pagein 
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algebraic (finite-difference) equations. Algorithms of this type have been widely used 
both in axisymmetry and in fully generic slices. 

3.2.5.1 Angular Coordinates, Grid, and Boundary Conditions 

In more detail, elliptic-PDE algorithms assume that the horizon is a Strahlkorper 
about some local coordinate origin, and choose an angular coordinate system and a 
finite-difference grid of points on in the manner discussed in section imi 

The most common choices are the usual polar-spherical coordinates {6, 0) and a 
uniform “latitude/longitude” grid in these coordinates. Since these coordinates are 
“unwrapped” relative to the actual trial-horizon-surface topology, the horizon 
shape function h satisfies periodic boundary conditions across the artificial grid 
boundary at 0 = 0 and 0 = 27r. The north and south poles 6 = 0 and 6 = tt are 
trickier, but Huq et al. [T^IHU]. Shibata and Uryu ina, and Schnetter j1 1 81 m h] all 
describe suitable “reaching across the pole” boundary conditions for these artificial 
grid boundaries. 

Alternatively, Thornburg |14()j avoids the z axis coordinate singularity of polar- 
spherical coordinates by using an “inflated-cube” system of 6 angular patches to 
cover S'^. Here each patch’s nominal grid is surrounded by a “ghost zone” of 
additional grid points where h is determined by interpolation from the neighboring 
patches. The interpatch interpolation thus serves to tie the patches together, en¬ 
forcing the continuity and differentiability of h across patch boundaries. Thornburg 
reports that this scheme works well but was quite complicated to program. 

Overall, the latitude/longitude grid seems to be the superior choice: it works 
well, is simple to program, and eases interoperation with other software. 

3.2.5.2 Evaluating the Expansion 0 

The next step in the algorithm is to evaluate the expansion 0 given by (EH) 
on the angular grid given a trial horizon surface shape function h on this same 
grid (jl.bj) . 

Most researchers compute 0 via 2-dimensional angular finite differencing of (El) 
on the trial horizon surface. 2nd order angular finite differencing is most common, 
but Thornburg jl4()j uses 4th order angular finite differencing for increased accuracy. 

With a (0,0) latitude/longitude grid the Q{h,duh,duvh) function in (13. 2 j) is 
singular on the 2 ; axis (at the north and south poles 0 = 0 and 0 = vr), but can 
be regularized by applying L’Hopital’s rule. Schnetter [TTR1ITTT1] observes that using 
a Cartesian basis for all tensors greatly aids in this regularization. 

Huq et al. [Zni EO] choose, instead, to use a completely different computation 
technique for 0, based on 3-dimensional Cartesian finite differencing: 

for further discussion of spectral methods.) This should yield a highly efficient apparent-horizon 
finder. However, I know of no published work taking this approach. 
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1. They observe that the scalar held F dehned by (EZD can be evaluated at any 
(3-diniensional) position in the slice by computing the corresponding (r, 6, 0) 
using the usual hat-space formulas, then interpolating h in the 2-dimensional 
{9, (f)) surface grid. 

2. Rewrite the apparent horizon condition (EH) in terms of F and it’s (3-dimensional) 
Cartesian derivatives, 

0 = 0(F, diF, dijF; gij, dkgij, Kij) = 0 (3.10) 

Huq et al. [ZgED] give the Q{F, diF, dijF) function explicitly. 

3. For each (latitude/longitude) grid point on the trial horizon surface, dehne a 
3x3x3-point local Cartesian grid centered at that point. The spacing of this 
grid should be such as to allow accurate hnite diherencing, i.e. in practice it 
should probably be roughly comparable to that of the underlying numerical- 
relativity simulation’s grid. 

4. Evaluate F on the local Cartesian grid as described in step Q above. 

5. Evaluate the Cartesian derivatives in dTTUH by centered 2nd order Cartesian 
hnite diherencing of the F values on the local Cartesian grid. 

Comparing the diherent ways of evaluating 0, 2-dimensional angular hnite dif¬ 
ferencing of (EH) seems to me to be both simpler (easier to program) and likely 
more efficient than 3-dimensional Cartesian hnite diherencing of (j3.1 Ojl . 

3.2.5.3 Solving the Nonlinear Elliptic PDE 

A variety of algorithms are possible for actually solving the nonlinear elliptic 
PDE (El (or (ITTUl) for the Huq et al. [HEni horizon hnder). 

The most common choice is to use some variant of Newton’s method. That is, 
hnite diherencing (El or (ITTUD (as appropriate) gives a system of A"ang nonlinear 
algebraic equations for the horizon shape function h at the Ng,ng angular grid points; 
these can be solved by Newton’s method in Aang dimensions. (As explained by 
Thornburg una section VIII.C], this is usually equivalent to applying the Newton- 
Kantorovich algorithm f |331 appendix C]) to the original nonlinear elliptic PDE (I3.2|l 
or (EH-) 

Newton’s method converges very quickly once the trial horizon surface is suffi¬ 
ciently close to a solution (a marginally outer trapped surface). However, for a less 
accurate initial guess, Newton’s method may converge very slowly or even fail to 
converge at all. There’s no usable way of determining a priori just how large the 
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radius of convergence of the iteration will be, but in practice ^ to | of the horizon 
radius is often a reasonable estimated^ 

Thornburg nsg described the use of various “line search” modifications to New¬ 
ton’s method to improve its radius and robustness of convergence, and reported 
that even fairly simple modihcations of this sort roughly doubled the radius of 
convergence. 

Schnetter [TTRIITTTI] used the PETSc general-purpose elliptic-solver library m 
EDIIZI] to solve the equations. This offers a wide variety of Newton-like algorithms 
already implemented in a highly optimized form. 

Rather than Newton’s method or one of its variants, Shibata et al. [TTmiTT^ use 
a functional iteration algorithm directly on the nonlinear elliptic PDE (Id.211 . This 
seems likely to be less efficient than Newton’s method but avoids having to compute 
and manipulate the Jacobian matrix. 


3.2.5.4 The Jacobian Matrix 


Newton’s method, and all its variants, require an explicit computation of the 
Jacobian matrix 

- II P'li) 

where the indices i and j label angular grid points on the horizon surface (or equiv¬ 
alently on S^). 

Notice that J includes contributions both from the direct dependence of 0 on h, 
duh, and d^vh, and also from the indirect dependence of 0 on h through the position- 
dependence of the geometry variables gij, dkQij, and Kij (since 0 depends on the 
geometry variables at the horizon surface position, and this position is determined 
by h). Thornburg jl37j discusses this indirect dependence in detail. 

There are two basic ways to compute the Jacobian matrix. 


Numerical Perturbation: The simplest way to determine the Jacobian matrix 
is by “numerical perturbation”, where for each horizon-surface grid point j, h is 
perturbed by some (small) amount e at the jth grid point (that is, hi ^ hi + eJu), 
and the expansion 0 is recomputed.The jth column of the Jacobian matrix (I3.11|l 
is then estimated as 


0i(h -f eJij) — 0i(h) 


'ij 


(3.12) 


^^Thornburg |137| used a Monte-Carlo survey of horizon-shape perturbations to quantify the 
radius of convergence of Newton’s method for apparent-horizon finding. He found that if strong 
high-spatial-frequency perturbations are present in the slice’s geometry then the radius of conver¬ 
gence may be very small. Fortunately, this problem rarely occurs in practice. 

very important optimization here is that 0 only needs to be recomputed within the finite 
difference domain of dependence of the Jth grid point. 
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Curtis and Reid m and Stoer and Bulirsch jl351 section 5.4.3] discuss the optimum 
choice of £ in this algorithm.^® 

This algorithm is easy to program but somewhat inefficient. It’s used by a 
number of researchers, including Schnetter j1 18| 11 19] and Huq et al. izniED]. 


Symbolic Differentiation: A more efficient, although somewhat more compli¬ 
cated, way to determine the Jacobian matrix is the “symbolic differentiation” algo¬ 
rithm described by Thornburg HSU, and also used by Pasch nng, Shibata et al. unD 
1132] , and Thornburg HU. Here the internal structure of the hnite differenced 0(h) 
function is used to directly determine the Jacobian matrix elements. 

This algorithm is best illustrated by an example which is simpler than the full 
apparent horizon equation: Suppose we discretize the left hand side L of the ap¬ 
parent horizon equation (ISZ3) with centered 2nd order hnite differences in 6 and 0. 
Then neglecting hnite-differencing trunation errors, and temporarily adopting the 
usual notation for 2-dimensional grid functions, hij = h{6=6i,(j)=(f)j), L is given by 


T _ hi—l,j 1 l,j 

^ (Ah)2 “^t^ 2Ah 


siffih 

The Jacobian of L is thus given by 

r 1 

± 


(A0)^ 


dL( 

dh, 


(Ah)2 2 tanh Ah 

1 


h’t) _ sin^ 9 (A0)2 


(fc/) 


if {k,i) = {i±l,j) 
if {k,i) = 
if {k,i) = (ij) 


(Ah)2 sin^h(A0)2 
0 otherwise 


(3.13) 


(3.14) 


Thornburg HSTI describes how to generalize this to nonlinear differential operators 
without having to explicitly manipulate the nonlinear hnite diherence equations. 


3.2.5.5 Solving the Linear Equations 

All the algorithms described in section in.2.5.3l for treating nonlinear elliptic PDEs 
require solving a sequence of linear systems of N^ng equations in N^ng unknowns. 
A^ang is typically on the order of a few thousand, and the Jacobian matrices in 

^®Because of the one-sided finite differencing, the approximation (Tna is only 0{e) accurate. 
However, in practice this doesn’t seriously impair the convergence of a horizon finder, and the 
extra cost of a centered-finite-differencing O(e^) approximation isn’t warranted. 
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question are sparse due to the locality of the angular hnite differencing (see sec¬ 
tion ld.2.5.4|l . Thus, for reasonable efficiency, it’s essential to use linear solvers 
that exploit this sparsity. Unfortunately, many such algorithms/codes only handle 
symmetric positive-definite matrices while, due to the angular boundary conditions^® 
(see section umi, the Jacobian matrices that arise in apparent-horizon finding 
are generally neither of these. 

The numerical solution of large sparse linear systems is a whole subfield of 
numerical analysis. See, for example. Duff, Erisman, and Reid ra and Saad |llbj 
for extensive discussions.^^ In practice, a numerical relativist is unlikely to write her 
own linear solver but, rather, will use an existing subroutine (library). 

Kershaw’s jSni ILUCG iterative solver is often used; this is only moderately effi¬ 
cient, but is quite easy to program.Schnetter [TTRIITTTI] reports good results with 
an ILU-preconditioned GMRES solver from the PETSc library. Thornburg |14()j 
experimented with both an ILUCG solver and a direct sparse LU decomposition 
solver (Davis’s UMFPACK library [22112211211120]), and found each to be more 
efficient in some situations; overall, he found the UMFPACK solver to be the best 
choice. 

3.2.5.6 Sample Results 

As an example of the results obtained with this type of apparent-horizon finder, 
figure IT21 shows the numerically-computed apparent horizons (actually, marginally 
outer trapped surfaces) at two times in a head-on binary black hole collision. (The 
physical system being simulated here is very similar to that simulated by Matzner 
et al. Eg, a view of whose event horizon is shown in figure |231) 

As another example, figure 13.31 shows the time dependence of the irreducible 
masses of apparent horizons found in a (spiraling) binary black hole collision, simu¬ 
lated at several different grid resolutions, as found by both AHFinderDirect and 
another CACTUS apparent-horizon finder, AHFlNDER.^® For this evolution, the two 
apparent-horizon finders give irreducible masses which agree to within about 2% for 
the individual horizons and 0.5% for the common horizon. 

As a final example, figure 12.71 shows the numerically-computed event and ap¬ 
parent horizons in the collapse of a rapidly rotating neutron star to a Kerr black 
hole. (The event horizons were computed using the EHFinder code described in 

the interpatch interpolation conditions in Thornburg’s multiple-grid-patch scheme USD]. 

^^Multigrid algorithms are also important here; these exploit the geometric structure of the 
underlying elliptic PDE. See Briggs, Henson, and McCormick m and Trottenberg, Oosterlee, and 
Schuller da for general introductions to multigrid algorithms. 

^^Madderom’s Fortran subroutine DILUCG m has been used by a number of numerical rela¬ 
tivists for both this and other purposes. 

^^AHFinder incorporates both a minimization algorithm fsection ft. 2.If II and a fast-flow algo¬ 
rithm lsection nf.2.7.411 : these tests used the fast-flow algorithm. 
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Figure 3.2: This figure shows the numerically-computed apparent horizons (actually 
marginally outer trapped surfaces) at two times in a head-on binary black hole 
collision. Figure reprinted with permission from Thornburg, Classical and Quantum 
Gravity 21, 743-766. Copyright 2004 by lOP Publishing Ltd. 
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dx=0.080 (AHFinder) 
dx=0.080 (AHFinderDirect) 


dx=0.060 (AHFinderDirect) 
dx=0.048 (AHFinderDirect) 
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Figure 3.3: This figure shows the irreducible masses (-^/area/lbvr) of individual 
and common apparent horizons in a binary black hole collision, as calculated by two 
different apparent-horizon finders in the Cactus toolkit, AHFinder and AHFind¬ 
erDirect. (AHFinderDirect was also run in simulations at several different 
resolutions.) Notice that when both apparent-horizon finders are run in the same 
simulation (resolution (ix=0.080), there are only small differences between their 
results. Figure reprinted with permission from Alcubierre et ai, Physical Review D 
72, 044004 (2005). Copyright 2005 by the American Physical Society. 
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section 


3.2.5.7 Summary of Elliptic-PDE Algorithms/Codes 

Elliptic-PDE apparent-horizon finders have been developed by many researchers, 
including Eardley EH, Cook jini UHl EZl , and Thornburg nsg in axisymmetry, and 
Shibata et al. irrmiTT^ . Huq et al. EH ED], Schnetter ITTRIITTTI] . and Thornburg |140j 
in fully generic slices. 

Elliptic-PDE algorithms are (or can be implemented to be) generally the fastest 
horizon-finding algorithms. For example, Thornburg jl4()j reports that the produc¬ 
tion version of his AHFinderDirect elliptic-PDE apparent-horizon finder, when 
run at each time step of a binary black hole evolution, averaged 1.7 seconds per 
time step, as compared with 61 seconds for an alternate “fast-fiow” apparent-horizon 
finder AHFiNDER (discussed in more detail in section rh2.7j) . However, achieving 
maximum performance comes at some cost in implementation effort (e.g. the “sym¬ 
bolic differentiation” Jacobian computation discussed in section in.2.5.4j) . 

Elliptic-PDE algorithms are probably somewhat more robust in their conver¬ 
gence (i.e. they have a slightly larger radius of convergence) than other types of 
local algorithms, particularly if the “line search” modiffcations of Newton’s method 
described by Thornburg nsTi are implemented.^"^ Their typical radius of convergence 
is on the order of 30% of the horizon radius, but cases are known where it’s much 
smaller. For example, Schnetter, Herrmann, and Pollney jl21j report that (with no 
“line search” modifications) it’s only about 10% for some slices in a binary black 
hole coalescence simulation. 

Schnetter’s TGRapparentHorizon 2D |1181 IllOj and Thornburg’s AHFind¬ 
erDirect jl40j are both elliptic-PDE apparent-horizon finders implemented as 
modules (“thorns”) in the Cactus computational toolkit. Both are freely available 
by anonvmous CVS. and work with either the PUGH unigrid driver or the Car¬ 
pet mesh-refinement driver for CACTUS TGRapparentHorizon2D is no longer 
maintained, but AHFinderDirect is actively supported and is now used by many 
different research groups.^® 

3.2.6 Horizon Pretracking 

Schnetter et al. mnmn] introduced the important concept of “horizon pretrack¬ 
ing”. They focus on the case where we want to find a common apparent horizon 
as soon as it appears in a binary black-hole (or neutron-star) simulation. While a 
global (flow) algorithm fsection Fd. 2. 7j) could be used to find this common apparent 

^^The convergence problems Thornburg ina noted when high-spatial-frequency perturbations 
are present in the slice’s geometry, seem to be rare in practice. 

^®In addition, at least two different research groups have now ported, or are in the process of 
porting, AHFinderDirect to their own (non- Cactus ) numerical relativity codes. 
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horizon, these algorithms tend to be very slow. They observe that the use of a local 
(elliptic-PDE) algorithm for this purpose is somewhat problematic: 

The common [apparent] horizon [... ] appears instantaneously at some 
late time and without a previous good guess for its location. In practice, 
an estimate of the surface location and shape can be put in by hand. 

The quality of this guess will determine the rate of convergence of the 
finder and, more seriously, also determines whether a horizon is found at 
all. Gauge effects in the strong field region can induce distortions that 
have a large influence on the shape of the common horizon, making them 
difficult to predict, particularly after a long evolution using dynamical 
coordinate conditions. As such, it can be a matter of some expensive trial 
and error to hnd the common apparent horizon at the earliest possible 
time. Further, if a common apparent horizon is not found, it is not 
clear whether this is because there is none, or whether there exists one 
which has only been missed due to unsuitable initial guesses - for a fast 
apparent horizon hnder, a good initial guess is crucial. 

Pretracking tries (usually successfully) to eliminate these difficulties by deter¬ 
mining - before it appears - approximately where (in space) and when (in time) the 
common apparent horizon will appear. 

3.2.6.1 Constant-Expansion Surfaces 

The basic idea of horizon pretracking is to consider surfaces of constant expansion 
(“CE surfaces”), i.e. smooth closed orientable 2-surfaces in a slice satisfying the 
condition 

Q = E (3.15) 

where the expansion E is a specihed real number. Each marginally outer trapped 
surface (including the apparent horizon) is thus a CE surface with expansion E = 0] 
more generally (IXT^ defines a 1-parameter family of 2-surfaces in the slice. As 
discussed by Schnetter et al. [TMITTg . for asymptotically flat slices containing a 
compact strong-held region, some of the E > 0 members of this family typically 
foliate the weak-held region. 

In the binary-coalescence context, for each t = constant slice we dehne E^ to be 
the smallest E > 0 for which a CE surface (containing both strong-held regions) 
exists with expansion E. If E* = 0 then this “minimum-expansion CE surface” is the 
common apparent horizon, while if E* > 0 this surface is an approximation to where 
the common apparent horizon will appear. We expect the minimum-expansion 
CE surface to change continuously during the evolution, and its expansion E^ to 
decrease towards 0. Essentially, horizon pretracking follows the time evolution of 
the minimum-expansion CE surface and uses it as an initial guess for (searching for) 
the common apparent horizon. 
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3.2.6.2 Generalized Constant-Expansion Surfaces 

Schnetter ing implemented an early form of horizon pretracking, which followed 
the evolution of the minimum-expansion constant-expansion surface, and found that 
it worked well for simple test problems. However, Schnetter et al. Qzn found that 
for more realistic binary-black-hole coalescence systems the algorithm needs to be 
extended: 

• While the expansion is zero for a common apparent horizon, it’s also zero 
for a 2-sphere at spatial inhnity. Figure 03] illustrates this for Schwarzschild 
spacetime. Notice that for small positive E* there will generally be two distinct 
CE surfaces with E = E*, an inner surface just outside the horizon, and an 
outer one far out in the weak-held region. The inner CE surface converges to 
the common apparent horizon as E* decreases towards 0, and is the surface we 
would like the pretracking algorithm to follow. Unfortunately, without mea¬ 
sures such as those described below, there’s nothing to prevent the algorithm 
from following the outer surface, which does not converge to the common 
apparent horizon as decreases towards 0. 
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Figure 3.4: This hgure shows the expansion 0 (left scale), and the “generalized 
expansions” r0 (left scale) and r^0 (right scale), for various r = constant surfaces 
in an Eddington-Finkelstein slice of Schwarzschild spacetime. Notice that all three 
functions have zeros at the horizon r = 2m, and that while 0 has a maximum at 
r 4.4m, both r0 and r^0 increase monotonically with r. 
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• In a realistic binary-coalescence simulation, the actual minimum-expansion CE 
surface may be highly distorted, which makes it hard to represent accurately 
with a finite-resolution angular grid. 

Schnetter et al. lEH discuss these problems in more detail, arguing that to solve 
them, the expansion 0 should be generalized to a “shape function” H given by one 
of 


Hi = Q 
Hr = he 
Hr^ = h‘^e 

CE surfaces are then generalized to surfaces satisfying 

H = E 


(3.16a) 

(3.16b) 

(3.16c) 


(3.17) 


for some specihed E > 0. 

Note that unlike Hi, both Hr and Hr^ are typically monotonic with radius. 
Neither Hr nor are 3-covariantly defined, but they both still have the property 
that E = 0 in ()3.17j) implies the surface is a marginally outer trapped surface, and 
in practice they work better for horizon pretracking. 

3.2.6.3 Goal Functions 

To select a single “smallest” surface at each time, Schnetter et al. jl21j introduce 
a second generalization, that of a “goal function” G, which maps surfaces to real 
numbers. The pretracking search then attempts, on each time slice, to hnd the 
surface (shape) satisfying H = E with the minimum value of G. They experimented 
with several different goal functions, 

Gh = H (3.18a) 

GrH — hH (3.18b) 

Gr = h (3.18c) 

where in each case the overline ( ) denotes an average over the surface.^® 

3.2.6.4 The Pretracking Search 

Schnetter’s uni original implementation of horizon pretracking (which followed 
the evolution of the minimum-expansion CE surface) used a binary search on the de¬ 
sired expansion E. Because E appears only on the right hand side of the generalized 

Schnetter et al. use a simple arithmetic mean over all surface grid points. In theory this 
average could be defined 3-covariantly by taking the induced metric on the surface into account, 
but in practice they found that this wasn’t worth the added complexity. 
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CE condition (ITT7I) . it’s trivial to modify any apparent-horizon finder to search for a 
snrface of specified expansion E. (Schnetter used his TGRapparentHorizon2D 
elliptic-PDE apparent-horizon finder described in section l!1.2.5.7l for this.) A binary 
search on E can then be used to hnd the minimum value 

Implementing a horizon-pretracking search on any of the generalized goal func¬ 
tions dSIHl) is conceptually similar but somewhat more involved: As described by 
Schnetter et al. inn for the case of an elliptic-PDE apparent-horizon hnder,^® we 
Erst write the equation defining a desired pretracking surface as 

H -H + G-p = 0 (3.19) 

where p is the desired value of the goal function G. Since H is the only term in (I3.19|l 
which varies over the surface, it must be constant for the equation to be satisfied. 
In this case H — H vanishes, so the equation just gives G = p, a.s desired. 

Because H depends on H at all surface points, directly hnite differencing 
would give a non-sparse Jacobian matrix, which would greatly slow the linear- 
solver phase of the elliptic-PDE apparent-horizon finder fsection l3.2.5.5j) . Schnetter 
et al. inn section III.B] show how this problem can be solved by introducing a 
single extra unknown into the discrete system. This gives a Jacobian which has a 
single non-sparse row and column, but is otherwise sparse, so the linear equations 
can still be solved efficiently. 

When doing the pretracking search, the cost of a single binary-search iteration is 
approximately the same as that of Ending an apparent horizon. Schnetter et al. w 
Egure 5] report that their pretracking implementation (a modiEed version of Thorn¬ 
burg’s AHFinderDirect ms elliptic-PDE apparent-horizon Ender described in 
section typically takes on the order of 5 to 10 binary-search iterations. 

The cost of pretracking is thus on the order of 5 to 10 times that of Ending a 
single apparent horizon. This is substantial, but not prohibitive, particularly if the 
pretracking algorithm isn’t run at every time step. 

^^There is one complication here: Any local apparent-horizon finding algorithm may fail if the 
initial guess isn’t good enough, even if the desired surface is actually present. The solution is to 
use the constant-expansion surface for a slightly larger expansion E as an initial guess, gradually 
“walking down” the value of E to find the minimum value i?*. Thornburg pni appendix C] 
describes such a “continuation-algorithm binary search” algorithm in detail. 

^®So far as I know this is the only case that has so far been considered for horizon pretracking. 
Extension to other types of apparent-horizon finders might be a fruitful area for further research. 

^®This refers to the period before a common apparent horizon is found. Once a common apparent 
horizon is found, then pretracking can be disabled as the apparent-horizon finder can easily “track” 
the apparent horizon’s motion from one time step to the next. 

^°With a binary search the number of iterations depends only weakly (logarithmically) on the 
pretracking algorithm’s accuracy tolerance. It might be possible to replace the binary search by 
a more sophisticated 1-dimensional search algorithm (I discuss such algorithms in appendix 0, 
potentially cutting the number of iterations substantially. This might be a fruitful area for further 
research. 
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3.2.6.5 Sample Results 

As an example of the results obtained from horizon pretracking, hsTure ITHl shows 
the expansion 0 for various pretracking surfaces (i.e. various choices for the shape 
function R in a head-on binary black hole collision. Notice how all three of the shape 
functions (IXTHD result in pretracking surfaces whose expansions converge smoothly 
to zero just when the apparent horizon appears (at about t = 1.1). 

As a further example, hgure ITHl shows the pretracking surfaces at various times 
in a spiraling binary black hole collision, projected into the black hole’s orbital plane. 

3.2.6.6 Summary of Horizon Pretracking 

Pretracking is a very valuable addition to the horizon-hnding repertoire: it es¬ 
sentially gives a local algorithm (in this case an elliptic-PDE algorithm) most of 
the robustness of a global algorithm in terms of hnding a common apparent horizon 
as soon as it appears. It’s implemented as a higher-level algorithm which uses a 
slightly-modihed elliptic-PDE apparent-horizon Ending algorithm as a “subroutine”. 

The one significant disadvantage of pretracking is its cost: each pretracking 
search typically takes 5 to 10 times as long as finding an apparent horizon. Further 
research to reduce the cost of pretracking would be useful. 

Schnetter et a/.’s pretracking implementation |121j is implemented as a set of 
modifications to Thornburg’s AHFinderDirect jl4()j apparent-horizon finder. Like 
the original AHFinderDirect, the modified version is a “thorn” in the Cactus 
toolkit and is freely available by anonymous CVS. 

3.2.7 Flow Algorithms 

Flow algorithms define a “flow” on 2-surfaces, i.e., they define an evolution of 
2-surfaces in some pseudo-time A, such that the apparent horizon is the A —cx) 
limit of a (any) suitable starting surface. Flow algorithms are different from other 
apparent-horizon finding algorithms (except for zero-finding in spherical symmetry), 
in that their convergence doesn’t depend on having a good initial guess. In other 
words, flow algorithms are global algorithms 1 section l3.1.7|l . 

To find the (an) apparent horizon, i.e., an outermost MOTS, the starting surface 
should be outside the largest possible MOTS in the slice. In practice, it generally 
suffices to start with a 2-sphere of areal radius substantially greater than 2mADM- 

The global convergence property requires that a flow algorithm always flow from 
a large starting surface into the apparent horizon. This means that the algorithm 
gains no particular benefit from already knowing the approximate position of the 
apparent horizon. In particular, flow algorithms are no faster when “tracking” the 
apparent horizon (repeatedly finding it at frequent intervals) in a numerical time 
evolution. (In contrast, in this situation a local apparent-horizon finding algorithm 
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Head-on collision 
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Figure 3.5: This figure shows the expansion 0 for various pretracking surfaces, 
i.e., for various choices for the shape function H, in a head-on binary black hole 
collision. Notice how the three shape functions (jd.lbjl (here labelled 0, r0, and 
r^0) result in pretracking surfaces whose expansions converge smoothly to zero just 
when the apparent horizon appears (at about t = 1.1). Notice also that these three 
expansions have all converged to each other somewhat before the common apparent 
horizon appears. Figure reprinted with permission from Schnetter, Herrmann, and 
Pollney, Physical Review D 71, 044033 (2005). Copyright 2005 by the American 
Physical Society. 
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Figure 3.6: This figure shows the pretracking surfaces at various times in a spiraling 
binary black hole collision, projected into the black holes’ orbital plane. Notice how, 
even well before the common apparent horizon first appears {t = 16.44mADM, bottom 
right plot), the r0 pretracking surface is already a reasonable approximation to the 
eventual common apparent-horizon’s shape. Figure reprinted with permission from 
Schnetter, Herrmann, and Pollney, Physical Review D 71, 044033 (2005), Copyright 
2005 by the American Physical Society. 
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can use the most recent previously-found apparent horizon as an initial guess, greatly 
speeding the algorithm’s convergence.) 

Flow algorithms were hrst proposed for apparent-horizon hnding by Tod |142j . 
He initially considered the case of a time-symmetric slice (one where Kij =0). In 
this case a marginally outer trapped surface (and thus an apparent horizon) is a 
surface of minimal area, and may be found by a “mean curvature flow” 

= -Ks^ (3.20) 

where x® are the spatial coordinates of a horizon-surface point, s® is as before the 
outward-pointing unit 3-vector normal to the surface, and k = is the mean 

curvature of the surface as embedded in the slice. This is a gradient flow for the 
surface area, and Grayson HH has proven that if the slice contains a minimum-area 
surface, this will in fact be the stable A —cx) limit of this flow. Unfortunately, this 
proof is valid only for the time-symmetric case. 

For non-time-symmetric slices. Tod |142j proposed generalizing the mean curva¬ 
ture flow to the “expansion flow” 


= -0s® (3.21) 

There is no theoretical proof that this flow will converge to the (an) apparent 
horizon, but since the flow velocity is zero there, and the flow is identical to the 
mean curvature flow in the principle part, convergence is at least theoretically 

plausible. Numerical experiments by Bernstein [21], Shoemaker et al. jl331ll34] . and 
others show that that the expansion flow does in fact converge robustly to 

the apparent horizon. 

In the following subsections I discuss a number of important implementation 
details for, and rehnements of, this basic algorithm. 

3.2.7.1 Implicit Pseudo-Time Stepping 

Assuming the Strahlkorper surface parameterization (D, the expansion flow 
is a parabolic equation for the horizon shape function This means that any 
fully explicit scheme to integrate it (in the pseudo-time A) must severely restrict its 
pseudo-time step AA for stability, and this restriction grows (quadratically) worse 
at higher spatial resolutions.^^ This makes the horizon-hnding process very slow. 

To avoid this restriction, practical implementations of flow algorithms use im¬ 
plicit pseudo-time integration schemes; these can have large pseudo-time steps and 
still be stable. Because we only care about the X —>■ oo limit, a highly accu¬ 
rate pseudo-time integration isn’t important; only the accuracy of approximating 

^^Linearizing the Q{h) function lit.211 gives a negative Laplacian in h as the principal part. 

^^For a spatial resolution Ax, an explicit scheme is generally limited to a pseudo-time step 
AA < (Aa:)h 
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the spatial derivatives matters. Bernstein m nsed a modified Dn Fort-Frankel 
scheme |SH] bnt fonnd some problems with the snrface shape gradually developing 
high-spatial-frequency noise. Pasch |l()8j reports that an “exponential” integrator 
(Hochbruck et al. [ZZj) works well, provided the flow’s Jacobian matrix is com¬ 
puted accurately. The most common choice is probably that of Shoemaker 
et al. jU.Jl 1184] . who use the iterated Crank-Nicholson (“ICN”) scheme.^® They 
report that this works very well; in particular, they don’t report any noise problems. 

By refining his finite-element grid (section ll.2.Jj) in a hierarchical manner, Met¬ 
zger [HHl is able to use standard conjugate-gradient elliptic solvers in a multigrid¬ 
like fashion,using each refinement level’s solution as an initial guess for the next 
higher refinement level’s iterative solution. This greatly speeds the flow integration: 
Metzger reports that the performance of the overall surface-finding algorithm is 
“of the same order of magnitude” as that of Thornburg’s AHFinderDirect jl4()j 
elliptic-PDE apparent-horizon finder (described in section irf.2.5.7|l . 

In a more general context than numerical relativity, Osher and Sethian have 
diseased a general class of numerical algorithms for integrating “fronts propagating 
with curvature-dependent speed”. These flow a level-set function 1 section 11.2.1|1 
which implicitly locates the actual “front”. 

3.2.7.2 Varying the Flow Speed 

Another important performance optimization of the standard expansion flow 
is to replace 0 in the right-hand side by a suitable nonlinear function of 0, chosen 
so the surface shrinks faster when it’s far from the apparent horizon. For example. 
Shoemaker et al. use the flow 


dxx^ 


(0 — c) arctan^ 


0 -c 


(3.22) 


for this purpose, where 0o is the value of 0 on the initial-guess surface, and c (which 
is gradually decreased towards 0 as the iteration proceeds) is a “goal” value for 0. 

^^Richtmyer and Morton ITTKl section 7.5] give a very clear presentation and analysis of the 
Du Fort-Frankel scheme. 

^^More precisely, Pasch |1 d8| found that that an exponential integrator worked well when 
the flow’s Jacobian matrix was computed exactly (using the symbolic-differentiation technique 
described in section However, when the Jacobian matrix was approximated using the 

numerical-perturbation technique described in section Pasch found that the pseudo-time 

integration became unstable at high numerical resolutions. 

•^^Pasch im also notes that the exponential integrator uses a very large amount of memory. 

•^^Teukolsky a and Leiler and Rezzolla [HJ have analyzed ICN’s stability under various 
conditions. 

■^^See the references cited in footnote El on page 1441 for general introductions to multigrid 
algorithms for elliptic PDFs. 
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3.2.7.3 Surface Representation and the Handling of Bifurcations 

Since a flow algorithm starts with (topologically) a single large 2-sphere, if there 
are multiple apparent horizons present the surface must change topology (bifurcate) 
at some point in the flow. Depending on how the surface is represented, this may 
be easy or difficult. 

Pasch jl()8j and Shoemaker et al. jl3311134j use a level-set function approach 
fsection ll.2.1j) . This automatically handles any topology or topology change. How¬ 
ever, it has the drawback of requiring the flow to be integrated throughout the 
entire volume of the slice (or at least in some neighborhood of each surface). This 
is likely to be much more expensive than only integrating the flow on the surface 
itself. Shoemaker et al. also generate an explicit Strahlkorper surface representation 
fsection ll. 2. 2|1 . monitoring the surface shape to detect an imminent bifurcation and 
reparameterizing the shape into 2 separate surfaces if a bifurcation happens. 

Metzger uses a finite-element surface representation fsection fl.2.311 . which 
can represent any topology. However, if the flow bifurcates, then to explicitly rep¬ 
resent each apparent horizon the code must detect that the surface self-intersects, 
which may be expensive. 

3.2.7.4 Gundlach’s “Fast Flow” 

Gundlach [72] introduced the important concept of a “fast flow”. He observed 
that the subtraction and inversion of the flat-space Laplacian in Nakamura et a/.’s 
spectral integral-iteration algorithm (section (21231) is an example of “a standard way 
of solving nonlinear elliptic problems numerically, namely subtracting a simple linear 
elliptic operator from the nonlinear one, inverting it by pseudo-spectral algorithms 
and iterating”. Gundlach then interpreted Nakamura et a/.’s algorithm as a type 
of flow algorithm where each pseudo-time step of the flow corresponds to a single 
functional-iteration step of the Nakamura et al. algorithm. 

In this framework, Gundlach defines a 2-parameter family of flows interpolating 
between Nakamura et a/.’s algorithm and Tod’s |142j expansion flow (EUD, 

dxh = -A{1- (3.23) 

where H > 0 and B > 0 are parameters, p > 0 is a weight functional which depends 
on h through at most 1st derivatives, is the flat-space Laplacian operator, and 
(1 — RV^)“^ denotes inverting the operator (1 — RV^).^® 

Gundlach then argues that intermediate “fast flow” members of this family 
should be a useful compromises between the speed of Nakamura et a/.’s algorithm, 
and the robustness of Tod’s expansion flow. Based on numerical experiments, Gund¬ 
lach suggests a particular choice for the weight functional p and the parameters A 

^®The inversion is only formal, because Nakamura et al.'s algorithm treats the ooo spectral 
coefficient specially. Gundlach m discusses this in more detail. 
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and B. The resulting algorithm updates high-spatial-frequency components of h 
essentially the same as Nakamura et a/.’s algorithm, but should reduce low-spatial- 
frequency error components faster. 

Alcubierre’s AHFinder |1] horizon hnder includes an implementation of Gund- 
lach’s fast flow algorithm.^® AHFinder is implemented as a module (“thorn”) in 
the Cactus computational toolkit, and is freely available by anonymous CVS (it’s 
part of the CactusEinstein set of thorns included with the standard Cactus 
distribution). AHFinder has been used by a large number of research groups. 

3.2.7.5 Summary of Flow Algorithms/Codes 

Flow algorithms are the only truely global apparent-horizon finding algorithms, 
and as such can be much more robust than local algorithms. In particular, flow algo¬ 
rithms can guarantee convergence to the outermost MOTS in a slice. Unfortunately, 
these convergence guarantees hold only for time-symmetric slices. 

In the forms which have strong convergence guarantees, flow algorithms tend 
to be very slow. (Metzger’s algorithm is a notable exception: it’s very fast.) 
There are modifications which can make flow algorithms much faster, but then 
their convergence is no longer guaranteed. In particular, practical experience has 
shown that in some binary black hole coalescence simulations (Alcubierre et al. [5], 
Diener et al. [Ml), “fast flow” algorithms fsection f3.2.7.4|l can miss common apparent 
horizons which are found by other (local) algorithms. 

Alcubierre’s apparent-horizon finder AHFinder ^ includes a “fast flow” algo¬ 
rithm based on the work of Gundlach |Z2|.“ It’s implemented as a module (“thorn”) 
in the Cactus computational toolkit, and is freely available by anonymous CVS 
(it’s part of the CactusEinstein set of thorns included with the standard Cactus 
distribution). It has been used by a number of research groups. 

3.3 Summary of Algorithms/Codes for Finding 
Apparent Horizons 

3.3.1 Summary of Apparent-Horizon Finding Algorithms 

There are a large number of apparent-horizon finding algorithms, with differing 
trade-offs between speed, robustness of convergence, accuracy, and ease of program¬ 
ming. 

In spherical symmetry, zero-finding (section 13.2.Ij) is fast, robust, and easy to 
program. In axisymmetry, shooting algorithms (section I3.2.2j) work well and are 

AHFinder also includes a minimization algorithm Isection 
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fairly easy to program. Alternatively, any of the algorithms for generic slices (snm- 
marized below) can be used with implementations tailored to the axisymmetry. 

Minimization algorithms (section Id. 2 .dll are fairly easy to program, but are sus¬ 
ceptible to spurious local minima, have relatively poor accuracy, and tend to be very 
slow unless axisymmetry is assumed. 

Nakamura et a/.’s spectral integral-iteration algorithm 1 section fd. 2 . 4jl and elliptic- 
PDE algorithms (section ld.2.5j] are both fast and accurate, but are moderately 
difficult to program. Their main disadvantage is the need for a fairly good initial 
guess for the horizon position/shape. 

In many cases Schnetter’s “pretracking” algorithm (section ld.2.tij) can greatly 
improve an elliptic-PDE algorithm’s robustness, by determining - before it appears 
~ approximately where (in space) and when (in time) a new outermost apparent 
horizon will appear. Pretracking is implemented as a modihcation of an existing 
elliptic-PDE algorithm, and is moderately slow: it typically has a cost 5 to 10 times 
that of Ending a single horizon with the elliptic-PDE algorithm. 

Finally, flow algorithms (section ld.2.7|l are generally quite slow (Metzger’s al¬ 
gorithm jHHl is a notable exception), but can be very robust in their convergence. 
They are moderately easy to program. Flow algorithms are global algorithms, in 
that their convergence does not depend on having a good initial guess. 

3.3.2 Summary of Publicly-Available Apparent-Horizon Find¬ 
ing Codes 

I know of d publicly-available apparent-horizon Ending codes, all implemented 
as modules (“thorns”) in the Cactus computational toolkit: 

AHFinder 

Alcubierre’s AHFinder |3] includes both a “fast flow” algorithm based on the 
work of Gundlach and a minimization algorithm based on the work of 
Anninos et al. [7j. AHFinder is part of the CactusEinstein set of thorns 
included with the standard Cactus distribution, and has been used by many 
research groups. 

AHFinderDirect 

Thornburg’s AHFinderDirect moi uses an elliptic-PDE algorithm, and 
has been used by many research groups, as well as ported to at least two non- 
Cactus numerical relativity codes. Schnetter’s pretracking algorithm inn is 
implemented as a set of modifications to AHFinderDirect. 

TGRapparentHorizon2D 

Schnetter’s TGRapparentHorizon2D jllHMllbj uses an elliptic-PDE algo¬ 
rithm. It’s no longer maintained, but remains freely available. 
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Appendix A 

Solving A Single Nonlinear 
Algebraic Equation 


In this appendix I briefly outline numerical algorithms and codes for solving 
a single 1-dimensional nonlinear algebraic equation f(x) = 0, where the function 
/ : 3? ^ 3ft is given. 

The process generally begins by evaluating / on a suitable grid of points and 
looking for sign changes. Assuming / to be continuous, each sign change must then 
bracket at least one root x^: 

Given a pair of ordinates a;_ and x^ which bracket a root, there are a variety of 
different algorithms available to accurately and efficiently find the (a) root: 

If |a;+ — X-\ is small, say on the order of a finite-difference grid spacing, then 
closed-form approximations are probably accurate enough: 

• The simplest approximation is a simple linear interpolation of / between a;_ 
and x+. 


• A slightly more sophisticated algorithm, “inverse quadratic interpolation”, is 
to use 3 ordinates, two of which bracket a root, and estimate the root as 
the root of the (unique) parabola which passes through the 3 given [x, f{x)) 
points.^ 

For larger \x+ — X-\, iterative algorithms are necessary to obtain an accurate 
root: 


• Bisection (binary search on the sign of /), is a well-known iterative scheme 
which is very robust, but rather slow if high accuracy is desired. 

• Newton’s method can be used, but it requires that the the derivative /' be 
available. Alternatively, the secant algorithm (similar to Newton’s method but 

^The parabola generically has two roots, but normally only one of them lies between x- and 
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estimating f from the most recent pair of function evaluations) gives similarly 
fast convergence without requiring /' to be available. Unfortunately, if |/'| is 
small enough at any iteration point, both these algorithms can fail to converge, 
or more generally they can generate “wild” trial ordinates. 

• Probably the most sophisticated algorithm is that of van Wijngaarden, Dekker, 
and Brent. This is a carefully engineered hybrid of the bisection, secant, and 
inverse quadratic interpolation algorithms, and generally combines the rapid 
convergence of the secant algorithm with the robustness of bisection. The 
van Wijngaarden-Dekker-Brent algorithm is described by Forsythe, Malcolm, 
and Moler |bdl chapter 7], Kahaner, Moler, and Nash |8dl chapter 7], and 
Press et al. mu section 9.3]. An excellent implementation of this, the Fortran 
subroutine ZEROIN, is freely available from http://www.netlib.org/fmm/ 
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Appendix B 

The Numerical Integration of 
Ordinary Differential Equations 


The time-integration problem^for ordinary differential 
traditionally written as follows: We are given an integer n 
ODEs to integrate), a “right-hand-side” function f : 3?” x 3? 
y(0) of a function y ; 3? —> 3fJ” satisfying the ODE 

| = /(y.O (B.i) 

We wish to know (or approximate) y(t) for some hnite interval t G [0,fmax]- 

This is a well-studied problem in numerical analysis. See Forsythe, Malcolm, 
and Moler j63l chapter 6] or Kahaner, Moler, and Nash |83l chapter 8] for a general 
overview of ODE integration algorithms and codes, or Shampine and Gordon [T2S1, 
Hindmarsh IZEI, or Brankin, Gladwell, and Shampine OH for detailed technical 
accounts. 

For our purposes, it suffices to note that highly accurate, efficient, and robust 
ODE-integration codes are widely available. In fact, there is a strong tradition 
in numerical analysis of free availability of such codes. Notably, the RKF45 
code described by Forsythe, Malcolm, and Moler ESI chapter 6] is freely available 



ODEPACK/LSODE suite of codes described by Hindmarsh [7^] are freely 

^The numerical-analysis literature usually refers to this as the “initial value problem”. Un¬ 
fortunately, in a relativity context this terminology often causes confusion with the “initial data 
problem” of solving the ADM constraint equations. I use the term “time-integration problem for 
ODEs” to (try to) avoid this confusion. 

^In this appendix sans-serif lower-case letters abc...z denote variables and functions in 3?" 
(for some fixed dimension n), and sans-serif upper-case letters ABC ... Z denote nx n real-valued 
matrices. 


equations (ODEs) is 
> 0 (the number of 
—>• 3fJ"', and the value 
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available at http://www.netlib.org/odepack/, and the RKSUITE suite of 
codes described by Brankin, Gladwell, and Shampine j211 are freely available at 
http://www.netlib.org/ode/rksuite/ 

As well as being of high numerical quality, these codes are also very easy to 
use, employing sophisticated adaptive algorithms to automatically adjust step size 
and/or the precise integration scheme used.^ These codes can generally be relied 
upon to produce accurate results both more efficiently and more easily than a hand¬ 
crafted integrator. I have used the LSODE solver in several research projects with 
excellent results. 


^LSODE can also automatically detect and handle stiff systems of ODEs. 
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